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Preface 


Perturbation theory is fun, useful, and, we believe, accessible to undergrad- 
uates in engineering and the physical sciences. The mathematical back- 
ground we expect of you, the reader, is modest: a two-semester course 
in one-variable calculus and a one-semester course in ordinary differential 
equations. Our book is written in an informal style, stressing heuristics. We 
have tried always to move from specific examples to generalities, emphasiz- 
ing the why along with the how.“ Both of us have used the material in 
this book in classes, and we know that the ideas can be grasped if you work 
consciously with them. This means that you should read this book with 
pencil and paper at hand to perform some of the implied algebraic opera- 
tions and that you must work most of the exercises. Perturbation theory, 
like any art, must be learned by doing. Fortunately, many of the tedious 
calculations in perturbation theory traditionally carried out by hand can 
now be relegated to the computer thanks to symbol manipulation programs, 
such as MACSYMA developed at the Laboratory for Computer Science at 
the Massachusetts Institute of Technology. 

Since perturbation methods produce only approximate solutions, one 
may ask, “Why not use numerical methods?” One answer is that per- 
turbation methods produce analytic approximations that often reveal the 
essential dependence of the exact solution on the parameters in a more sat- 
isfactory way than does a numerical solution. A second answer is that some 
problems which cannot be easily solved numerically may yield to pertur- 
bation methods. Indeed, numerical and perturbation methods should be 
combined in a complementary way. 

Let us mention the organization of the book. Chapter I is an overture 
that introduces the major themes that are elaborated upon in the chap- 
ters that follow. Chapter II is a thorough, but not exhaustive, treatment 
of how to find the roots of polynomials whose coefficients contain a small 
parameter. This chapter introduces regular and singular perturbations in 
as simple a context as possible. Working through this chapter carefully 
will help you to fix the idea of perturbation methods while using only al- 
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gebraic manipulations to solve problems. Subsequent chapters concentrate 
on differential equations. Here, we introduce you to many techniques for 
handling perturbations which change the order of the equations or which 
work for differential equations whose independent variable is unbounded. 
We end the book with two disparate practical problems that can be solved 
efficiently with perturbation methods. 

Many people have influenced this book, but it is with special warmth 
that we recall the lessons of two of our respective teachers-George Carrier 
(jem) and Eric Reissner (jgs)-who introduced us to the power and delights 
of perturbation theory. And we remember the deep and dazzling mastery 
of the subject by our friend and former colleague, Gordon Latta. 

Finally, our thanks to Carolyn Duprey whose efforts with the computer 
made multiple revisions of the manuscript possible. 


Since writing the first edition of this book, we have become more con- 
vinced not only for the need for an elementary introduction to this material 
but also for the need to examine all mathematical formulations of models 
of natural phenomenon to see if small or large parameters allow the struc- 
ture of the formulation to be simplified. Indeed, the odds are overwhelming 
that if any physical system is nondimensionalized the resulting dimension- 
less parameters will be very small or very large, so it makes sense to develop 
methods to exploit this windfall. As computing power increases, it becomes 
easier to “get an answer” without understanding important aspects of the 
governing laws. The thinking required to understand perturbation theory 
stands against the presentation of numerical computation without thought 
for the meaning or correctness of such calculations. If critical examina- 
tion of results can be learned as an undergraduate, much frustrating and 
needless effort can be saved in the workplace and in graduate school. 

In this revision, we have sought to correct errors and to provide slightly 
different wording in some places. We have also recognized that in a course 
on perturbation theory the instructor might like to raise the topic of ap- 
proximation of integrals. To this end, we have included an appendix to 
introduce this important topic. 


J. G. S. 
Charlottesville, Va. 
J. E. M. 

Wheaton, II. 


Chapter 1 


Introduction and 
Overview 


Perturbation theory is the study of the effects of small disturbances. If 
the effects are small, the disturbances or perturbations are said to be regular; 
otherwise, they are said to be singular. The basic idea in perturbation 
theory is to obtain an approximate solution of a mathematical problem by 
exploiting the presence of asmall dimensionless parameter—the smaller the 
parameter, the more accurate the approximate solution. 

Regular perturbations are assumed nearly every time we construct a 
mathematical model of a real world phenomenon. Our choice of language 
reflects this: the flow is almost steady, the density varies essentially with 
altitude only, the conductivity is virtually independent of temperature, the 
spring is nearly linear, the friction is practically negligible. 

Singular perturbations are probably less familiar. Fig. 1.1 illustrates 
two examples. The top in Fig. I. Ia is set spinning rapidly about a vertical 
axis. During one revolution, the effects of aerodynamic drag and the friction 
at the tip are small (regular perturbations). Eventually, though, the top 
falls and comes to rest in a position far from its initial one. Thus, over a 
long period of time, the perturbations are singular. Models such as this are 
characterized by what may be called a singularity in the domain. For the 
top, this means that we are interested in what happens for all time after its 
release, that is, for all times t on the semi-infinite (and therefore singular) 
domain t > 0. 

Fig. 1.1b shows a hemispherical elastic shell under an internal pressure 
p. The shell has been clamped at its edge, which prevents displacement or 
rotation there. Since the shell is symmetric, the stresses depend only on 
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b) 


Figure 1.1: Singular perturbations resulting from (a) a singularity in the 
domain and (b) a singularity in the model. 


the polar angle ¢. The maximum stress occurs in the outer fibers of the 
shell and can be found by solving non-homogeneous ordinary differential 
equations. The solutions contain constants determined by regularity con- 
ditions at the pole (¢ = 0) and edge conditions at the equator (¢ = 7/2). 
The differential equations contain the small parameter h/R, where h is the 
shell thickness and R is the radius of the shell midsurface. If h/R is set to 
zero in these equations, then we obtain the equation for a membrane , i.e., 
a shell with no bending stiffness. The solutions of these simplified equa- 
tions predict a maximum stress of pR/2h everywhere.! The defect of these 
simplified equations is that their solutions cannot meet the edge conditions. 


This result is also derived easily if we divide a complete spherical membrane by an 
imaginary equatorial plane and then equate the pressure force, prR?, to the tension 
along the equator, 2r Rho, where ø is the stress. 
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This fact signals a singularity in the model. 

In Fig. 1.1b we have plotted, as a function of ¢, the maximum dimen- 
sionless stress (2ha/pR) as predicted by shell and membrane theory for a 
typical value of h/ R. Except near the equator, the results are virtually iden- 
tical. However, in a narrow zone near the equator—the boundary layer—the 
dimensionless maximum stress predicted by shell theory dips below 1 and 
then rises to a value of 2. The key feature of this graph is that no matter 
how small the parameter h/ N, the rise of the stress by a factor of 2 at the 
edge never diminishes. This is a singular perturbation phenomenon. The 
“shell versus membrane” solutions reflect what may be called a singularity 
in the model. Setting h/R = 0 leads to an over-simplified model that fails 
to predict the non-negligible stress rise at the boundary. This failure of the 
membrane theory occurs in a narrow region near the boundary; the width 
of the failure region depends on the size of h/R. 

Perturbation problems arising from a singularity in the domain were 
first studied systematically by Poincaré, who encountered them in orbital 
mechanics. The first extensive analysis of problems involving a singular- 
ity in the model (boundary layer problems) was done by Prandtl in his 
study of low viscosity fluids flowing over solid objects. Although the prob- 
lems attacked by Poincaré and Prandtl are too elaborate for this book, we 
can explore many aspects of perturbation theory by working with simple 
equations, many of which can be applied to common phenomena. 

The First Quantitative Step. We begin with the following problem 
from the theory of quadratic equations: 


Determine how the roots of z? — 2z + change as € is perturbed 
slightly away from zero. 


If e = O, the roots are, by inspection, 21 = 0 and 22 = 2. For other 
values of e we have, as a result of the quadratic formula, 


a(e)=1—-Vl-e (1.1) 


za(e)=1+Vl-e. (1.2) 


The Numerics of a Regular Problem. Using a hand calculator, ? we 
readily construct from (1.1) and (1.2) the following table. Our numerical 
calculations suggest that a perturbation about e€ = 0 is regular but teach 
us little else. Moreover, we took no advantage of the fact that the roots 
for e = O came with little effort. We shall remedy this lack of analysis 
presently. 


2A Texas Instruments TI-30-II in our case 
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Table 1.1. Roots of z? — 2z +€ 


z2(e) 
1.9486833 - - - 


1 19949874 - 
10-6 | 0.000005 --- | 1.999999. 
10-7 | 5.01 x 1078 | 1.9999999... 


The Numerics of a Singular Problem. If we switch e with the 
coefficient of 22, we have the polynomial ez? — 2z +1. Ife = 0, z1 = 1/2 is 
its only root. If e # 0, there are two: 


20 — (1.3) 


ad= 2 — = (1.4) 


Again, using a hand calculator, we — the following table. What is 


Table 1.2. Roots of ez? — 2z ＋ 1 


€ zı(€) z2(€) 
10-! | 0.5131670--- | 19.4868333 - -- 
107? | 0.5012563---| 199.49874-.- 


100.500 1000 .. 1999999.5--- 
10 7 0.501 19999999 


notable about Table 1.2? First, while one of the roots, 21 (e), approaches 
that of 2z — 1 as e 0, the other, 22 (6), goes to infinity. This is a manifes- 
tation of singular behavior. On a more fundamental level, changing e from 
0 to an arbitrarily small number has changed the number of solutions for 2 
(from one to two). Notice that the equation changes from linear (e = 0) to 
quadratic (e # 0). Such a change in the order of an equation characterizes 
many singular perturbation problems. 

A second feature of Table 1.2 is that as € — 0, 21 (e) at first approaches 
0.5 steadily but then, for very small values of €, begins to exhibit small 
fluctuations. These are caused by round-off errors produced by computing 


A First Look At Perturbation Theory 5 


differences of nearly equal terms in (1.3). While there is a simple way to 
remedy this in the present problem—multiply numerator and denominator 
in (1.3) by 1+ /1 — e—the diagnosis of round-off error in more complicated 
problems, much less the cure, is not so simple. In cases such as these in 
which numerics falter, perturbation theory can sometimes save the day. 

Analysis of the Regular Problem. To find approximate solutions 
that are accurate and easy to use, we study the effect of a small parameter. 
First, let us find approximate formulas for the roots of z? — 2z + e when 
€ is small. Unlike most perturbation problems, this one can be analyzed 
completely. We analyze precisely the simplest problem of a class with the 
hope of inferring patterns or principles which can aid in the attack on more 
complicated problems. 

Infinite Series. The formulas (1.1) and (1.2) for the roots of 22 22 
are exact. However, the smallness of e simplifies the problem. We can 
certainly assume that fe < 1. This not only rules out complex roots, but, 
more importantly, it allows us to expand /1—€ in a power series in e. 
Recall the binomial expansion: 


= ! 
(a+b)™ = am ma . * HY Fr — 
(1.5) 


Setting a = 1, b e, and m= 1/2, we have 


WI Ee=I- -= (1.6) 


This is a formal series, so-called because we as yet have made no attempt 
to ask what it means to try to add together an infinite number of terms. 

The way to study an infinite series is to study its sequence of partial 
sums. If the sequence converges, then we say that the series converges, and 
if we are able to compute the limit, S, of this sequence, we say that the 
series sums to S. One of the simplest tests for convergence is the ratio test, 
which’ states that Lu converges if 


lim l <1. 
* 0 UE 
From (1.5) and (1.6) we find that 
„ |(e+1)*term| _ 1421/2 
E (1.7) 


Thus the right side of (1.6) converges if |e] < 1, as do the two series 


1 1 
2106) 1 1 +e 
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20 2 30 (1.8) 


that we obtain upon substituting (1.6) into (1.1) and (1.2). 

The first few terms of these series must, to be useful, yield close ap- 
proximations to the roots. Taking only those terms displayed explicitly 
in (1.8) and e = 0.1, we find 21 (.1) œ 0.5(.1) + 0.125(.1)? = 0.05125 and 
220.1) % 2 — 0.5(.1) = 1.95. The exact values are z;(.1) = 0.05131... and 
z2(.1) = 1.94.... Are these roots accurate enough? This decision is made 
from external information in the context of your own problem. 

In a more elaborate problem, there may not be a formula for the exact 
solution(s) with which one can check the accuracy of an approximation, or, 
if there is an exact formula, it may be difficult to apply. As an example of 
the latter, consider finding the roots of 


z3 — 522 T 42 T 2 0 (1.9) 


when ¢ is small. Though there is an exact formula for the roots of this (and 
any other) cubic, it is difficult to use. However, if € = 0, we have 


2 — 522 ＋ 42 = 0 (1.10) 
which has the roots 
2100) =0 
220) S 1 (1.11) 
230) = 4. 


Presently, we shall develop a technique which exploits the smallness of € to 
produce an acceptable approximation of the roots of any polynomial P,(z) 
when the roots of the nearby“ or reduced polynomial Po(z) are known. 
But first we need a way of making error estimates that depends on the 
approximation process itself. 

Taylor’s expansion with a remainder uses the value of a function 
and its first few derivatives at a point at which they are easily found to 
estimate the value of the function at a nearby point where the function 
is difficult to calculate. More precisely, if a function f(e) and its first n 
derivatives are continuous on a closed interval |e| < a, and if the (n+1)st 
derivative of f(e) exists on the open interval |e| < a, then 


(n) 
fle) = lo) e+ EQ ¹ο , (1.12) 
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Here F) (o) denotes the nth derivative of f(e) evaluated at € = 0 and the 
remainder after n+1 terms, Rn+1(€), has the form 


+e 
Rays) = FD at, al < Ie. (1.13) 


The number z depends on € but is otherwise unknown. 
The important facts about the expansion (1.12) are: 


1. it allows us to approximate f(€) by a (Taylor) polynomial in €, namely 
the right side of (1.12) without Nai (e). 


2. if we can find an upper bound on R,4:(e) then we have an upper 
bound on the error of our polynomial approximation. 


3. it is almost incidental in applications whether the associated infinite 
series f(0) + f’(0)e+ 1/2f"(0)e? +--- converges. The size of Rn41(€) 
is the salient fact. 


With f(e) = (1 — «)!/?, we have, ifn = 1, 
— r)-3/2 
(120% =1- x — Go, |z| < lel <1 (1.14) 


and, if n= 2, 


_ 7)-5/2 
(1-6/2 1 3- 5? — 1 — 4. |z| < le] <1. (1.15) 


Substituting (1.15) into (1.1) and (1.14) into (1.2), we find that 
a(e) = x + 5? + Cai lz| < ļeļ|< 1 


n- |z| < Jel <1. (1.16) 


These expressions offer a way of estimating the errors we make when trun- 
cating the infinite series expansions for z1 (e) and z2(e€) after two terms. 

If p > O, and |z| < 1, then (1 — z)~? is largest when z is as close to 
1 as possible. Thus, if |z| < 1 — ô, where ô is any fixed number such that 
0<6<1, 


(12) —½ , 

. 1857 leg. (1.17) 
123 37/2 

0 ) e? < 5 lel’. 
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Exercise 1.1. Earlier, we used the terms displayed explicitly 
in (1.8) to approximate z;(.1) and z2(.1). Use (1.16) and (1.17) 
to obtain an upper bound on the errors we made. 


The Order Symbols. Using (1.17), we may rewrite (1.16) in the form 


z(e) = x + 8 + Oles) (1.18) 
22 (60 = 2 — xt O(e?). (1.19) 


The symbols O(es) and O(e?), to be read “big ‘O’ of e cubed” and “big ‘O’ 
of e squared” are used to sweep all irrelevant algebraic details under the 
rug. 

In general, g(€) = O(e?) means that, for e sufficiently small, there exists 
a positive constant K independent of e, such that ge) < KJe|?. In (1.18) 
and (1.19), “sufficiently small” means |e| < 1—6, and the K’s from (1.17) are 
(1/16)5 5/2 and (1/8)6~/? respectively. In more complicated problems, 
however, we can rarely pin down the words “sufficiently small” and “there 
exists.” Thus, in practice, we may have to view a statement such as f(€) = 
O(e?) as simply implying that f grows no faster than the square of € when 
€ is small. 

Analysis of the Singular Problem. With what we have learned we 
can now quickly analyze the singular problem of finding simple, approxi 
mate formulas for the roots of ez? — 2z + 1 when e is small. Substituting 
the Taylor formula for VI € into (1.3) and (1.4), we obtain 


216) ; + e + Ole?) (1.20) 
oe z z ; +0(6). (1.21) 


The terms displayed give an approximation to 21010) of 0.5 + .125(107°) 
= 0.5000000125, and an approximation to z2(10~°) of 2(10°) — 0.5 = 
1999999.5. The approximation to 21 (10-5) has been improved dramatically 
from the value found earlier with a calculator, but 22 (e) still approaches 
infinity as € approaches zero. This behavior is inherent in the problem and 
is not a numerical artifact. Our analysis of 22(e) has, nevertheless, provided 
useful information: we now know that z2(e) behaves like 2/e as € — 0. For 
larger values of € we might need more terms in the Taylor polynomial for 
21(€) and z2(€) to obtain sufficiently accurate approximations. 
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Exercise 1.2. Make upper bound estimates of the remain- 
ders and determine the smallest Taylor polynomials that can 
produce approximations to the roots of ez? — 2z + 1 with an 
absolute error < 1073 for all |e| < 0.2. 


Note that when we finally obtained useful numerical formulas for the 
roots of z? —2z +e (1.18) and (1.19) each was of the form 


z(e)= ao+aje Tae f aN ENT RAI (e), Ryyi(e) = Ore et) . (1.22) 


The right side of (1.22) is called an asymptotic or regular perturbation ex- 
pansion. It is ideal for assessing, numerically or theoretically, the effect 
of a small perturbation in e about zero. Though an asymptotic expansion 
need not converge as N — oo, we do require that eN Ry = O ase 70 
for fixed N. Any function with a representation of the form (1.22) is called 
regular because z(€) approaches a finite value as € approaches 0. 

Re-analysis of the Regular Problem. Suppose we pretend that the 
quadratic formula does not exist. Let us assume instead that each root of 
P.(z) = z? — 2z + € has a representation of the form (1.22) for some fixed 
integer N. Then we may attempt to determine the unknown coefficients 
ao, 1 by requiring that 


P,(z(e)) = P(ao + aye 4A. + aye + O(e%+1)) 
= (ao tae +--+)? 
—2(ao ＋T die..) Te = 0, (1.23) 


identically in € as e — 0. In this problem we know that N may be any posi- 
tive integer, that ao +a,€+ dz? +- is a convergent power series, and that 
its radius of convergence is 1. However, as we shall see, such information 
is not necessary to determine the unknown coefficients in (1.22). Indeed, 
there are problems involving a small parameter in which the solutions are 
of the form (1.22), but in which the associated infinite series ao + det 
does not converge for any non-zero value of e. [For example, see (1.74)] 
Therefore, it is essential to emphasize that the subsequent procedure does 
not involve infinite series. 

Before substituting the right side of (1.22) into (1.23), we must square 
it. Thus, if z(€) = ao + Ri(e), where Ri(€) = Ole), then z?(e) = aĝ + 
2ao Ri (6) T Ri (e). Now 2aoR1(€) = O(e), because Ri (e) = O(e) implies that 
there exists a constant K such that |Ri(e)| < Ke| for e sufficiently small. 
Hence, 2d RI (e) < 24H el is true for |e| sufficiently small. Furthermore, 
Ri(e) = O(e?), because Ri (e)] < Ke] implies that |R2(e)| < K?le|?. In 
fact, since K?|e|? < K?|e| if |e] < 1, we can also write Ri (e) = O(e). Finally, 
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it is important to note that the sum of two O(e)-terms is again O(e). For 
if there exist constants P and Q such that |p(e)| < Ple| and |q(e)| < Ole 
for |e| sufficiently small, then |p(e) + ge) < Kle|, where K = P+Q. In 
summary, we may conclude that if z(€) = do +O(e), then 2206) = a3 Oe). 
See what a marvelous invention the O-symbols are! 

Again, if z(e) = ao + aieR2(€), where Ra(e) = O(e?), then 


27(€) = ag + Jade + aĵe? + 240 Na(e) + 2416 Ra(e) + Ne) 
= a3 + 2apa,e + O(e?). (1.24) 


To understand how we could make such a simplifying leap, let us take, 
one by one, the terms in the first line of (1.24) that we swept into the 
O-symbol in the second line. If K is any constant greater than |a?|, then 
|aze?| < K ſeſz, i.e., aĵe? = O(e?). (To allow for the possibility that the 
unknown coefficients might be complex numbers, we write |a?| instead of 
af.) The term 20 Ra(e) is also O(e?) because |R2(€)| < K|e|? implies that 
|2a0R2(€)| < 2|ao|K|e|?. For the same reason, 2a,¢R2(e) = O(e*). But a 
term that is O(es) is also O(c?) because Ke < Kle|? if |e] < 1. Finally, 
as |R2(€)| < K |e|? for e sufficiently small, it follows that 20 < K?\e|* < 
K?\e|? for |e| sufficiently small. That is, R3(e) is both O(e*) and O(e?). 
Thus the last four terms in the first row of (1.24) are each O(e?). But, by 
the same argument we used earlier to show that the sum of two O(e) terms 
is again O(e), we conclude that a finite sum of O(e?) terms is again O(c”), 
and thus arrive at the second line of (1.24). 

Before proceeding further into perturbation theory, be sure that you 
have a clear understanding of how O-symbols work; they are essential to 
the rest of this book. 

After computing explicitly one more term in the representation for 22 (e), 
we can write 


z7(€) = a} + 2apaye + (a? + 2apaz)e? + O(ê). (1.25) 
Substituting this expression and (1.22), with N = 2 into P,(z) = z?—2z+e, 
we obtain 
P.(z(€)) = a3 + 200g + (af + 20e + OC) 
—2[ao + aye + ane? Oe =0. (1.26) 
Collecting coefficients of like powers of € and noting that the sum of two 
terms of O(e*) is again a term of O(¢*), we can rewrite (1.26) in the form 
P,(z(€)) = aĝ — 2ao + (2aoa; — 2a; + 1)e 
+a? + 2000 — 2a2)€? + O(e*) = 0. (1.27) 
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It is important to keep in mind that to say that (1.22) is a root of P,(z), 
for e sufficiently small, means that P.(z(e)) must be identically zero (= 0) 
for all values of e sufficiently small. In particular, then, (1.27) must hold 
as € — 0. Therefore, because the left side is a continuous function of e, 


a2 — 2ao = 0, (1.28) 
This result in turn, reduces (1.27) to 
P.(z(e)) = (2aoa; — 2a; + 1) e + (a? + 2aoaz — 2az)e? + O(c) = 0. (1.29) 


But if P.(z(e)) = 0 for e sufficiently small, then e 1 P.(z(e)) = 0 for e 
sufficiently small and nonzero. This implies that lim e~1P,(z(e)) = 0 (“lim” 
is shorthand for lim. o) which, by (1.29), yields 


2041 — 2a, +1 = 0. (1.30) 
Thus (1.29) reduces still further to 
P.(z(e)) = (a? + 2aoa2 — 2a2)e? + O(e*) =0. (1.31) 
Again, since (1.31) must hold for all values of é sufficiently small, we con- 
clude that lime 2 B. (z(e)) = 0. This requires that 
az + 24002 — 202 = 0. (1.32) 


What we have just established, in a specific case, is 
The Fundamental Theorem of Perturbation Theory: 
If 


Ao + Are . . + Ave + O(e%*") = 0 (1.33) 


for e sufficiently small and if the coefficients Ap, A1, are inde- 
pendent of e, then 


A0 = Aj =- -= An =0. (1.34) 
To solve z? — 2z + e = 0, note that (1.28) has two roots that may be 
found by inspection: 


ao = 0 and do = 2. (1.35) 
Choose the first and (1.30) yields 
41 = 1/2. (1.36) 


Thus, from (1.32), 
(1/2)? — 2a = 0 or az = 1/8. (1.37) 


12 Introduction and Overview 


Substituting these results into (1.22), we reproduce (1.18), the representa- 
tion for 21 (€) obtained earlier from the quadratic formula. 

If we make the second choice which satisfies (1.28), ao = 2, and follow 
the same steps, we find that a1 = —1/2 and az = —1/8. These results, sub- 
stituted into (1.22), reproduce (1.19), the earlier representation for z2(e). 
The same procedures can be followed to produce as many terms of the se- 
ries for 21 (€) and 22 (e) as desired. 


Exercise 1.3. Use the above procedure to compute 21 (e) and 
2z2(€) when N = 4 in (1.22). 


Merely knowing that a solution to a problem exists is usually of little 
help in finding it. However, an important principle is illustrated by our 
example. If we can learn something about the general form of a solution, 
say that it has a representation like (1.22), then the unknown parts of the 
solution may often be determined one at a time by direct substitution into 


the given equations. 


Exercise 1.4. The roots of the polynomial z? — z + € can be 
shown to be of the form (1.22). Find ao and a; for each root. 


Re-analysis of the Singular Problem. Which roots of Q,(z) = 
ez? — 2z +1 have the form (1.22) with N = 1? Substituting (1.22) along 
with (1.25) for z?(e) into Q,(z), we find that 


Q.(z(€)) = ela + 2acaie + O(e?)] — 2[ag + aye + O(e?)J} +1. (1438) 


By collecting like powers of € and requiring that Q,(z(€)) be identically zero 
for e sufficiently small, we obtain 


- 2ao + 1 + (ag - 2a; )e + O(c?) =0. (1.39) 


Using The Fundamental Theorem of Perturbation Theory, we set the coef- 
ficients of the various powers of € in (1.39) to zero to obtain 


—2a9+1=0 
ag — 2a; = 0. (1.40) 
Solved in order, these equations yield 
40 = 1/2 
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a, = 1/8. (1.41) 


Substituting these values back into (1.22), we reproduce (1.20), the expres- 
sion we obtained earlier using the quadratic formula. What about the other 
root, z2(€), given by (1.21)? Obviously, we did not pick it up because (1.21) 
is not of the form (1.22). 

Recall that the point of our re-analysis is to pretend that we know 
neither the quadratic formula nor how to expand I- e in an infinite 
series. With this in mind, we can still infer that the missing root z2(e), 
must approach infinity as € approaches zero. Why? Because 


Q.(z) = ez? — 22 71 
= dz- (I= 2200 (1.42) 
= ez? — e[zi(€) + z2(€)]z + ezı (€)z2(€). 


Equating the 1 in the first line of (1.42) to €z;z2 in the last line of (1.42) 
(since the two lines must be equal if z = 0), we see that 


21 (c 22(6) = +. (1.43) 


This implies that as € — 0, zz (e) — oo, since z;(€), a regular root, must 
approach a finite value. 
Now we can make a key simplification; when z2(e) is large, —2z2(€) + 
1 ~ —2z2(e). [By f(e) ~ g(e), which is read “f is asymptotic to g as € 
approaches zero”, we mean that /g — 1 as € .] Thus 
Q.(z2(€)) ~ €23 (e) — 2z0(€) = He(z2(6)). (1.44) 


By inspection, the roots of H,(z) are 0 and 2/e. But since z(e) gets large 
as € gets small, it is evident that 


zale) ~ 2/e. (1.45) 


This result, obtained by heuristics (suggestive but not rigorous argu- 
ments), leads to the conjecture that 


czzl(e) ~ bo 4 0. (1.46) 


At this point, the fact that bọ = 2 is not essential; that bọ # 0 is. This 
conjecture in turn suggests the change of variable 


ez(e) = w(e). (1.47) 
Noting that eO. (z) = €?z? — 2ez + €, we have 
€Q.(z) = «Q(w/e) = w? — 2w + € = Ww). (1.48) 
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By (1.47) and (1.48), the roots of We (w) give the roots of Q,(z). But as 
Wo(z) has two roots, finding the roots of W,(w) for e sufficiently small is a 
regular perturbation problem. 


The essence of singular perturbation theory is to extract (by 
heuristics, if necessary) the dominant singular behavior of a 
solution and then, by a change of variable, to reduce the singular 
problem to a regular one. 


For the case at hand, we assume that each root of W,(w) has a regular 
perturbation expansion of the form (1.22) with, say, N = 1: 


(e) = bo +b16 + O(¢?). (1.49) 


Substituting (1.49) into (1.48), collecting coefficients of like powers of € and 
requiring that the resulting expression be identically zero for e sufficiently 
small, we find 


(b3 — 250) + (250b1 — 2b; + 1)e + O(c?) ~ 0. (1.50) 

By the fundamental theorem of perturbation theory, 
55 — 250 = 0 (1.51a) 
25001 — 2b; +1=0. (1.51b) 


The conjecture (1.46) implies, for the root we seek, that w(0) # 0. Thus, 
the solution of (1.51a) we want is 


50 = 2. (1.52a) 


From (1.51b) 
1 (1.52) 


The root of Re(z) we seek therefore has the representation 
oe x +0(€), (1.53) 
which, as a result of (1.47), produces 
2 1 
z2(e) = am 2 ＋ Ole). (1.54) 


Thus, we have reproduced (1.21) without using the quadratic formula. 
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Exercise 1.5. Find the first three non-zero terms in the ex- 
pansions of each of the roots of €z? — z? +1. 


The analysis of these two simple examples has suggested several basic 
principles, including the idea of reducing singular problems to regular ones 
by a change of variable and then seeking solutions of the form (1.22). But, 
are the solutions of every regular problem regular? The answer is, “not 
necessarily,” as the following problem shows. 


Find an expansion for the roots of P,(z) = z? — 2ez e, useful 
if |e| is small. 


The problem seems regular since no roots are lost if € = 0. Therefore, 
we seek solutions of the form (1.22). Following the familiar routine, we 
substitute (1.22) into P,(z), collect coefficients of like powers of e, and 
require that the resulting expression be identically zero if € is sufficiently 
small. With, say, N = 2, this yields 


a3 + (2aoa; — 2ao — 1)e 


+(a? + 240 — 241) e: + O(e*) = 0. (1.55) 
By the Fundamental Theorem of Perturbation Theory , 
ag =0 (1.56a) 
2aoa; — 2a9 —1 = 0 (1.56b) 
a? + 200% — 2a; = 0. (1.56c) 


The solution of (1.56a) is ao = 0 so that (1.56b) reduces to —1 = 0. This 
contradiction tells us there is no root of 22 — 2ez — € of the form (1.22). 

We could fall back on the quadratic formula, but the purpose here is 
to develop a line of reasoning that does not appeal to formulas for exact 
solutions. (Why? Because, for arbitrary polynomials of degree higher than 
four, no such formulas exist.) Instead, we start from the fact that the roots 
of 

P,(z) = 22 — 2ez — € (1.57) 

must approach zero as € does. We next ask how fast does this process 
occur? Let us assume that as € — 0, 


z(e) ~ ebo, p> 0, bo #0. (1.58) 


The symbol ~ means here that lim e Pz (e) = bo. Of course, there are many 
other ways that z(e) might approach zero, say z(e) ~ elne, but first we 
want to consider a simple power dependence. 
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Exponent 


Figure 1.2: Graphs of the exponents of ¢ in (1.60). Exponents are equal at 
‘intersections. 


The assumption (1.58) suggests the change of variable 


z(e) = ue), w(0) £0. (1.59) 
Substituting (1.59) into (1.57), we have 
P. (er w(e)) = «7? w?(e) 26e) — e. (1.60) 


As e 0, the largest of the three terms e, er, el is the one with the 
smallest exponent. If p > 1/2, then, as is clear from Fig. 1.2, min{2p, p+ 
1,1} =1. Hence, e~! P,(e?w(e)) ~ —1. But this is a contradiction since the 
left side must be identically zero for e sufficiently small. Likewise, we get 
a contradiction if p < 1/2 because in this situation min{2p, p+ 1,1} = 2p 
which implies that e 27 P.(€? w(e)) ~ w?(0) # 0 The only possibility left is 
p= 1/2. In this case, 
E! P(e Pw(e)) = w (e) — 2e/? w(e) — 1 ~ u - 1, (1.61) 
which implies that w(0) = +1. 
Observe in (1.61) that it is «/? that appears and not é. This suggests 


that, in addition to the change of variable (1.59), we introduce the change 
of parameter 


p=eP. (1.62) 


Thus, the problem of finding the roots of z? — 2ez —e for e sufficiently small 
has been replaced by the equivalent problem of finding the roots of 


So(w) = w? — 2Bw-1 (1.63) 


for @ sufficiently small. If our analysis has been correct, then Sg (u) should 
have two regular roots of the form 


w() = bo +618 T. +b 8% (N). (1.64) 
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Substitution of (1.64) into (1.63) leads to the sequence of equations 


53 — 12 0 
25061 — 250 = 0 
b? + 2bob2 — 251 =0 (1.65) 
whose solutions are, sequentially, 
bo = +1 
51 =1 
ba = +1/2. (1.66) 


Substituting into (1.64), (1.62), and (1.59), we find that the two roots of 
our original polynomial z? — 2ez c are, with N = 2 in (1.64), 


ale) / e+ 2 + O(e?) (1.67) 


zale) 402 +6- se? + O(¢?). (1.68) 


Exercise 1.6. Verify these results by using the quadratic for- 
mula [and show that the error terms are actually smaller than 
indicated in (1.67) and (1.68)}. 

Exercise 1.7. Determine explicitly the first two terms in the 
expansions of the three roots of 2° — ¢z? — 62. 


Exercise 1.8. Consider the linear algebraic equations 


er 2c T Z C 
—z—e*y+ez=1 
excteytz=e. 


As € = 0, z ae, y ~ doe’, z ~ coe”. Assume regular expan- 
sions of the form e = ao die , etc., and use Cramer’s 
Rule to determine p, q, r, ao, bo, co, a1, 61, c 

Exercise 1.9. Determine the dominant behavior as € — 0 of 
the eigenvalues of the coefficient matrix of the system of equa- 
tions in Exercise 1.8. 
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Accuracy of a Regular Expansion. Suppose that we can identify the 
first N +1 terms in (1.22) with the first N +1 terms of a power series. If the 
series has a radius of convergence p # 0, then, if |e| < p, we can compute 
z(e) to any degree of precision by making N large enough. Indeed, the 
examples considered so far have been of this type with p = 1. 

It may happen, however, that p = 0, in which case z (e) can be computed 
from (1.22) only to within some irreducible error. That is, for each value 
of e, it may happen that the remainder term in (1.22) can be made only 
so small. If we compute terms beyond a certain point, the approximation 
becomes less and less accurate. How many terms should be kept? With 
no information about the remainder term except that it is O(e"+1), the 
rule of thumb is: stop computing terms when the (n+1)st term in (1.22) is 
larger than the nth term. 

The following is a classic example of the situation described above, but 
may be skipped with no loss of continuity. 

The Small € Expansion of Stieltjes’ Integral 


«+f eae 
lle) = ji i (1.69) 


An expansion for € < 1 should not be too hard to find as J(0) = Sy etie 
1. To account for the effect of € in (1.69), € must be put into the numerator 


of the integrand. 
The following formula for a geometric series with a remainder is of great 
use. 
— =1l+z2+274+-- * + — . (1.70) 


Exercise 1.10. Verify (1.70). 


Setting z = —e*t and substituting (1.70) into (1.69), we have 


I(e) = 15 1 . (—1)"e?"4"] dt 


H-P PE iE), (1.71) 
where eee 
TIn+1(€) -f EI (1.72) 


Each of the integrals in the first line of (1.71) can be evaluated by techniques 
of calculus: 


co 
{ e™ttPdt =n! (1.73) 
0 
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Thus 
Ie) = 1 - e + 21e — . . 4 (1) ne 4 (1) 162 f() 
= S(-1)tetet# + Rnti(e). (1.74) 
0 


This expansion is of the same form as (1.22) with e2 in place of € and N 
= 2n. To show that it is regular, we need only bound Ini (e) by a constant 
independent of e. But 1 +e°t > 1 ift > 0. Hence, 


© tnt oo r i 
Bað [ Sear Sj T d= GY Ui, QD. (1.75) 


To show that Co (1) Ele cannot approximate I(e) to arbitrary pre- 
cision, no matter how we choose n, we use the last link in the following 
chain of lower bounds on In 1 (e): 


e dt 
1 40 EET (1.76) 


eit 
z f t + et 
co 
178 0 e (14 s) nds, t=1l+s 
0 


er! 7 
-8 N 
> 17 c sd 


With this result and (1.75), we obtain from (1.74), 
enn 


e(1+ €?) 
< IK) — }0(-1)* ke] (1.77) 


< enn + 1)! = Un41(€)- 


Laa (e) = 


Since 

Un(e) en = [e?(n)][e?(n — 1) ---[e*(2)]fe*(1)], (1.78) 
e?"n! begins to grow once ezn— the largest factor on the right side of 
(1.78)—exceeds one. Therefore, given e, en will be smallest when n 
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is the largest integer n. such that e°n, < 1. Moreover, as the right side 
of (1.78) is the product of n decreasing factors, its value is greater than 
the nth power of the smallest factor: that is, min en > (e) n > 2 if 
|e] < 1. It now follows from the left side of (1.78) that 


c2(t+e7?) 
Ln+1(€) > e+e) (1.79) 


Thus, we make an error at least as large as the right side of (1.79) when we 
approximate I(e) by Oo (1) E Ele. 

In practice, we would also like to know how small we can make Unie). 
It follows from (1.78) that, given e, Unie) is minimized by taking n as the 
largest integer such that n+1 < e~?. For example if e 5, taken = 3. This 
guarantees that I(.5) is approximated by 1—(.5)?+2(.5)*—6(.5)® = 0.78125 
to within an absolute error that is no greater than (.5)°4! = 0.09375. (The 
correct value of I(.5) is 4e4#1(4) = 0.82538 ., where EI (2) = f° e~*/t dt 
is the exponential integral. Values for Ei may be found in Abramowitz and 
Stegun, pp. 245-251.) 


Exercise 1.11. Graph the upper bound on the error in the 
optimum approximation to I(e) for 0 < € < 1. 
Exercise 1.12. Given e use the double inequality? 


Varn” t1/2 exp[—n + (12n + 1) 1 (1.80) 
<n! 
< Vian" +l? exp[—n (12n) 11 


to obtain lower and upper bounds on Ly41(€) and Uni (e), re- 
spectively, in terms of € only. According to these estimates, 
if e = .1, then )7$°(—1)*&!(.01)* should approximate I(.1) to 
within an error of 0(10~*°). 


Asymptotic Sequences. As we have tried to show by simple ex- 
amples, the goal of perturbation theory is to reduce singular problems to 
regular ones. These we try to solve approximately by solving a sequence of 
simpler problems that have been purged of all small parameters. Working 
different problems has forced us to keep enlarging our concept of a regular 


3See Feller, An Introduction to Probability and Its Applications, Vol. I, 3rd Ed., Wiley, 
1968, p. 54. 
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expansion (solution). At first, a regular expansion was merely a power se- 
ries in e with a known radius of convergence. Then, it was a power series 
convergent for some value of e sufficiently small. Next we decided that a 
finite power series with a remainder estimate might be a better definition 
of a regular expansion. When we saw that an equation containing integral 
powers of € could have a solution involving fractional powers of e, we again 
broadened our definition. 

Our final definition is this: a regular expansion is an expression of the 
form 


z(e) = doe) +41 Ge) +... + Gngn(e)+Rnsi(6), Rngile) = aie 


where 
lim $n41(€)/¢n(€) = 0 as € 0. (1.82) 


Any sequence of functions having the property (1.82) is said to be an asymp- 
totic sequence. You should convince yourself that the Fundamental The- 
orem of Perturbation , Theory, stated in the box containing (1.33) and 
(1.34), continues to hold if && is everywhere replaced by S (e). 

For the simple problems considered in the rest of the book, expansions 
of the form (1.22) are sufficient, but it is important to note that expansions 
of the form (1.81) are not uncommon. For example suppose that we need 
the small € expansion of 


fl Set elk. (1.83) 
From the power series expansion for the exponential function, we have 
€ 1 2 
e 1e + (1.84) 


(1e — ceeln € 


=I telne+ zen e+). (1.85) 


Thus 
f(E) = 14 % +e? Ine + (1/2)e? + (1/2)E In? + --- 
= Hole) + 261 (€) + d2(e) + (/) (e) + (1/2) (e) . (1.86) 
Clearly, the ¢,,’s satisfy (1.82). 


Differential Equations. Though finding the roots of a polynomial or 
the value of an integral is necessary in the analysis of a variety of problems, 
it is rare that a polynomial or integral itself models a phenomenon. (We 
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often try to fit experimental data with polynomials, but this is not true 
modelling. Rather, it is an admission that the mechanism producing the 
data is either too obscure or too complicated to have predictive value.) 
More often, differential equations (DE’s) are used as models, and for many 
phenomena, such as the motion of a point-mass in a gravitational field, 
the kinetics of a chemical reaction, or the growth of a colony of bacteria, 
ordinary DE’s (ODE’s) suffice. Our main concern shall be to determine 
when and how such ODE’s can be solved approximately by perturbation 
methods. 

Linear, constant coefficient ODE’s are encountered everywhere. In the- 
ory, they are easy to solve. It is important to remember that if the inde- 
pendent variable is, say t, then the first step toward a solution is to look 
for solutions of the form e, where z is a constant to be determined. In the 
case of a single equation of degree n (or a system of n first-order equations), 
this leads to finding the roots of a polynomial in z of degree n. If the coef- 
ficients of the DE (or DE’s) contain a parameter e, so will the polynomial, 
and if e is small, we may construct close approximations to the roots by 
using some of the techniques mentioned above. These techniques will be 
developed systematically in Chapter 2. 


Special Features of ODE’s. The solutions of ODE’s are subject to 
auxiliary conditions, usually in the form of initial conditions or boundary 
conditions, which may depend on the parameter €. More important, if the 
roots of the associated polynomial of a constant coefficient ODE depend on 
e, then the solution will depend on two variables, e and t. The interplay 
of these is crucial. A small disturbance, as reflected by the presence of a 
parameter e in a DE or its auxiliary conditions, may lead to a large effect 
because, when the range of t is unbounded, terms such as et can grow 
large (singularity in the domain) , or, if the range of t is bounded, say 
0 St < 1, terms such as t/e can grow large (singularity in the model). 
These phenomena are examples of nonuniformities. 


Nondimensionalization of the equations describing a physical system 
is done by expressing each variable having physical dimensions as a fraction 
of some fixed intrinsic quantity. Nondimensionalization is an extremely 
useful procedure, even when a problem is not amenable to perturbation 
methods. 

From physics we know that the period of a simple, frictionless pendu- 
lum of length /, oscillating with an arbitrarily small angular amplitude, is 
2n(I/9)'/?, where g is the strength of the gravitational field. Thus, in study- 
ing large oscillations, it is natural to measure the actual time as a fraction 
of (1/g)'/?, especially if we are interested in motions in which the amplitude 
and friction are small. The fraction itself is called the dimensionless time. 
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As another example, the equation of motion of a tightly-stretched string 
of length L takes a particularly concise form if distance along the string is 
represented by zL and time by tL(p/T)'/2. Here p is the mass per unit 
length, assumed constant, and T is the tension. Physically, L(p/T)¥ 2 is 
the time it takes a small disturbance to move from one end of the string to 
the other. 

Often there may be more than one way to nondimensionalize the equa- 
tions of a mathematical model. Consider the equations describing the fall 
of a ball released from rest at a distance h above the ground. Assuming 
a constant gravitational field of strength g and an atmosphere of constant 
density p, we could express the velocity as a fraction of (2gh)'/?, the veloc- 
ity of impact in vacuum. On the other hand, we could express the velocity 
as vV, where V, determined by equating the drag on the ball to the grav- 
itational force, is the terminal velocity the ball would approach if released 
from an infinite height. If the influence of the drag on the impact velocity 
is expected to be small, it is reasonable to choose (2gh)!/? as the unit of 
velocity; otherwise V is perhaps the best choice. 

A great virtue of nondimensionalization is that it is possible to represent 
an infinity of physical problems by a single mathematical problem. For 
example, the equation for the angular motion, (t), of a simple frictionless 
pendulum of any length / in a gravitational field g, released from rest at an 
angle €, is 4 - 

sin = 0, t>0, 600) Se, 0(0)=0. (1.87) 


Here, t is the dimensionless time discussed earlier and 6 = d6/dt. 
Exercise 1.13. Derive the DE in (1.87). 


Because of conservation of energy, the pendulum will return to its initial 
angular position after some dimensionless period of time t., which depends 
on e. The period of oscillation of the physical pendulum is (I/g)'/7t,. Thus, 
for a given initial angular displacement e, computation of the single number 
t. gives the period of a physical pendulum of any length / swinging in a 
gravitational field of any strength g.* 


Actually, from the principles of dimensionless analysis alone, it may be concluded 
that the period of a simple frictionless pendulum must be of the form (1/g)!/2t,(e). The 
function t. (e), however, can be determined only by solving (1.87), which we shall do in 
Chapter IV. A more striking example of the power of dimensional analysis comes from 
subsonic aerodynamics where it is shown that, when viscosity is negligible, the drag on 
an object immersed in a uniform flow of infinite extent is given by CHL Vn, where p 
is the air density, L is some characteristic dimension of the object, and V is the velocity 
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Even if effects such as friction are included, a single solution of the di- 
mensionless equation of motion still represents an infinite number of phys- 
ical solutions. For instance, suppose that our pendulum oscillates in a 
vacuum, but that its pivot exerts a torque 79 opposing the motion. Here ņ 
is a constant with the dimensions of (mass) x (length)? /(time). In place 
of (1.87), we have 


64+ 26 ＋ sin g = 0, t>0, 600) e, 6(0)=0, (1.88) 


where 2a = n/(mg!/?/°/?) and m is the mass of the bob of the pendulum. 
The period of motion is no longer constant. However, for given values of e 
and a, (1.88) represents an infinite class of equivalent physical problems 
two physical problems being equivalent if they have equal values of q and e. 


Exercise 1.14. An elastic beam of section modulus EI, resting 
on an elastic foundation of modulus k, is under a tension T and 
a distributed downward force/length p(s), where s is distance 
along the beam measured from some convenient point. The 
small vertical deflection w of the beam satisfies the DE 

d? dw dw 

ia ( 14%) T4 tius. (1.89) 
For simplicity, assume that EI, T, and k are constants, ex- 
pressed in some common set of physical units. Show, by setting 
s = Lz and appropriately choosing L (which has dimensions of 
length), that (1.89) can be given the dimensionless form 


(“ö Y y" =), (1.90) 
where y' = dy/dz. 


Another useful feature of nondimensionalization is that it sometimes 
reveals disparate physical systems to be identical mathematically. (This 
opens the possibility of solving physical problems by analogy.) Consider 
the forced, damped, linear spring-mass system in Fig. 1.3. Its motion is 
described by the differential equation and initial conditions 

2 
mos + 695 +k€=F(r), O<r (1.91) 


of the flow far from the object. The dimensionless drag coefficient Cp depends only 
on the shape of the object. Cp is extremely difficult to compute from the equations 
of fluid dynamics, even for the simplest shapes, and is therefore usually determined 
experimentally from wind tunnel tests. See Rauscher, M., Introduction to Aeronautical 
Dynamics, Wiley, 1953, Sections 10.7 & 10.8. 
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Figure 1.3: A linear, damped, spring-mass system. 


(%) Eo, €(0) = 0. 


(Other initial conditions could have been chosen). Here € is the location of 
the mass measured from the rest position of the spring, 7 is the time, m is 
the mass, c is a damping coefficient, and k the spring constant. The units 
of (e, r, m, c, k} are, say, {centimeters, seconds, grams, gram per second, 
gram per second squared}. 

The study of differential equations reveals that homogeneous solutions 
of the DE in (1.91) in the absence of damping (c = 0) are of the form 
cos(T\/k/m),sin(r,/k/m). Thus, the period of the free oscillations of the 
linear spring-mass system is 27,/m/k. Suppose that we are interested in 
weak damping forces. Because these will little change the undamped period 
of the system, it is natural to introduce the dimensionless time 


t TVA / m. (1.92) 


Furthermore, damping will cause the system to dissipate energy continually. 
In view of the initial conditions, the maximum displacement can never 
exceed o in magnitude if there is no forcing term. This suggests that we 
introduce the dimensionless displacement 


z = £/£o. (1.93) 


(Different initial conditions would suggest a different way of nondimension- 
aling €). With the change of variables (1.92) and (1.93), (1.91) is trans- 
formed into 

ž+2et +z = f(t), 0<t (1.94) 
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Figure 1.4: A linear LRC circuit. 


z(0)=1, 40) = O, 


where 


26 / VMR. (1.95) 


Now consider the closed circuit sketched in Fig. 1.4. that consists of a 
linear inductor, a linear resistor, and a linear capacitor in series and driven 
by an applied, time-varying voltage V(T)-a so-called LRC circuit. The 
current J in the circuit satisfies the differential equation 


dI dI I dv 
LTT RJ G 0<T (1.96) 


and a set of initial conditions, which we take to be 
1(0) = Io, IT) = 0. 


Here L, R, and C are, respectively, the inductance, resistance, and the 
capacitance of the system, measured, say, in units of Henrys, Ohms, and 
Farads, respectively. 

We can show that with the change of variables 


t=T/VLC, 2 = I/Io (1.97) 
and the introduction of the dimensionless parameter and function 
Sun, pal, (1.98) 
Io dT 


(1.96) reduces precisely to (1.94). That is, the forced motion of a damped 
linear oscillator is identical, mathematically, to the current in a driven LRC 
circuit. Find the solution for one of these systems and the solution of the 
other comes free. 
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Figure 1.5: An initially slack membrane under a central load. 


Exercise 1.15. As indicated in Fig. 1.5, the outer edge of an 
annular membrane is attached to a rigid plate with a hole of 
radius b. The inner edge of the membrane is attached to a mov- 
able rigid disk of radius a which is under a downward vertical 
force P. The membrane is of constant thickness h and obeys a 
linear, isotropic stress-strain law. According to Föppl's theory, 
the equilibrium and compatibility conditions for the membrane 
(assuming that fibers undergo small but finite angles of rota- 
tion) can be reduced to the single equation 


ala 1 EhP? 


where T is the radial tension and E is Young’s modulus, an 
elastic parameter here taken constant. 


(a) Show with the change of variables 


rz, f=¢T, (1.100) 
that (1.99) takes the form 
d? 1 EhP? 
4 -i r (1.101) 
(b) Set 
e S dir, FSA (1.102) 


and determine the constant À so that (1.101) takes the 
dimensionless form 
1 
* * (1.103) 
(e) State Newton’s Law for a mass-point falling towards 
the center of attraction of an inverse square gravita- 
tional force. Nondimensionalize the equation so that 


it reduces to (1.103). 


Chapter 2 


Roots of Polynomials 


In Chapter 1 we looked at several simple polynomials with coefficients de- 
pending on a parameter e. In each case we found that constructing a useful 
expansion for the roots reduced, ultimately, to computing, one-by-one, the 
coefficients in a regular expansion of the form (1.22). In some cases, all of 
the roots were regular. In other cases, a change of variable and perhaps 
parameter was necessary to get an associated polynomial that did have 


regular roots. 
The material of this chapter is not essential for the sequel, but here 


we develop a procedure for representing the roots of any polynomial of the 
form 


P.(z) = (1 + boe + coe? . ) + Are™(1 + bie. 92 (2.1) 

+--+ Ane (1 T bne r ), 
where the a;’s are rational numbers and the b;’s, c;’s, --- are constants. 
Further, each of the factors 1+ be T.., k = 0,1,--+,n, is assumed to 


be regular, i.e., to have an expansion of the form (1.22). The polynomials 
mentioned in Chapter 1 are special cases of (2.1). 


Exercise 2.1. Replace each of the following polynomials by a 
polynomial of the form (2.1) having the same non-zero roots: 
(a) e- 22 +2? 
(b) 2+6z+e7124 
(c) ez- 28. 
As e — 0, the different roots of P.(z) may approach zero, finite non- 
zero values, or infinity. For example, 1 + e~1z + 122 + 28 has roots of 


29 


30 Roots of Polynomials 


each type, as you’ll be asked to show later in an exercise. In Chapter 1 
we studied the roots of ez? — 2z + 1 and z? — 2ez — e. We found that the 
first polynomial had a non-regular root z2(€) ~ 261 and that the second 
had two non-regular roots, 21 (e) / and z2(e) ~ —e1/?. Expansions for 
these non-regular roots were obtained by introducing the change of variable 
z(e) = €?w(e) and then determining p by requiring that w(0) # 0. This is 
the key to finding useful expansions for the roots of (2.1) and motivates the 


following 
THEOREM: Each root of (2.1) is of the form 


z(e) = ), w(0) £0, (2.2) 


where w(e) is a continuous function of e for e sufficiently small. 

PROOF: First, we shall determine an algorithm for finding what we 
shall call the proper values of p. Then we shall establish (2.2) and the 
continuity of w(e). 

If (2.2) is to hold, then, from (2.1), 


P. (er u) = Oe u) + €(bo +t ́,E . eb Anw”)+---, (2.3) 


where 
Qe(w) = 1+ eM FP Aw H --- tP Anu". (2.4) 


Let z(e) = ue) be a root of P,(z), i.e., Pe(e w(e)) ~ 0 for all e. As 
Q.(w) is a continuous function of w and also of € if € # 0, it follows that if 
we) is continuous, then 


lim Q.(w(0)) = 0. (2.5) 


In view of (2.4), (2.5) implies that if w(0) # 0 then at least two of the 
exponents in the set 


E = {0, a1 +P, „ Qn + np} (2.6) 


must have identical, minimal values. The reason is clear: If w(0) # 0, 
then as € — 0 the dominant terms in Q,(w) are those whose ¢-factors 
have the smallest exponents. If there were but one minimal exponent, 
say a, + kp, then, since P,(ew(e)) = 0, it would follow from (2.3) that 
lim e~(¢*+#?) . (e w(e)) = Awt (0) = 0. But this is a contradiction, be- 
cause if a, + kp belongs to the set (2.6), Ay #0. 

To select the proper values of p, look at Fig. 2.1, the linear graphs of 
the exponents in (2.6). The circled intersections correspond to values of p 
where two or more exponents have equal, minimal values. If p is sufficiently 
large, the smallest exponent in the set (2.6) will be 0. As p decreases 
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1% 55% 


Pem) 


Figure 2.1: Graphs ofthe exponents of e in equation (2.4). 


(imagine a vertical line in Fig. 2.1 moving from right to left), there will 
be a first, largest proper value p; such that at least two exponents in (2.6) 
have the value ei = 0. One (and only one) of the graphs passing through 
the point (pi, ei) will have a largest slope nı. Now continue to decrease 
p until a second proper value pz is reached for which at least two of the 
exponents in (2.6) take on some identical, minimal value e2. The slopes of 
the graphs through (pz, ez) range from a minimum, ni, to some maximum, 
n2. Proceed in this fashion until the last and smallest proper value, pm, is 
reached. For p = pm at least two exponents in (2.6) will have the common 
value em = Q@n+Npm. Of the graphs through (pm, em), one has a minimum 
slope, nm-1, and that of a, + np has a maximum slope, nm = n. 

To make sure that it is clear how to compute the proper values of p, let 
us interrupt the proof of the theorem with an example, 


Pilz) = 1+ 23 + 626 + 2627 + 6228 + 6829 (2.7) 

(rigged, of course, to give simple answers.) The substitution z = eu yields 

PAP w) = 1+ Pw? + OP s 4 269. 7⁰ 77 4 c12 465 8 4 cls ere, (2.8) 
so the set of exponents of e is 

E = {0, 3p, 6+6p, 9+7p, 12+ 8p, 18 + 9p}. (2.9) 


From the graphs displayed in Fig. 2.2 we read off the proper values p; 
and the corresponding minimal exponents e;. These pairs are listed below 


along with the associated forms TË Xw) = e™ĉi P,(e?iw), whose importance 
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-35 
Figure 2.2: Graphs of the exponents of e in equation (2.8). 


will become clear in a moment. For the four choices of p, (O, —2, —3, —6), 
the scaled polynomials are: 


(0,0) : T) (w) = €°P.(e°w) = 14+ w te w2 w e ew? (2.10) 
(—2,-6) : TO (w) Se P.(e- 2) = w? + w® + 2ew + ew? (1 + u’) 


(2.11) 
(—3, —12) : T&)(w) = P ew) = w? + Qw7 + w? + (u? +w?) + e? 
(2.12) 
(—6, —36) : TP (w) = è P(w) (2.13) 
= u? +w? + 2607 + ew? + 618% + 686. 
Exercise 2.2. Given a set of rationals, [O, , an], write a 


program (using your favorite language or pseudo language) to 
compute the associated proper values of p. 


Back to the proof of the theorem. With the aid of (2.2), we rewrite (2.1) 
in the form 


T) = “P (Piw) = TP (w) EY), (2.14) 


where 


TË (w) = A, b +--+ Byw"s-*) (2.15) 


and EG (u) = 0. That is, Tw) has no e's and E%)(w) is zero when e 
is zero. Equations (2.10)—(2.14) are arranged in this form. What we have 
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done in multiplying P, by e~* and changing variables from z to w is to 
display explicitly the dominant part of P, as a polynomial TË ) which is 
independent of e. Clearly, TË ) has nj n -x non-zero roots. 

In Appendix A we prove the following plausible 

LEMMA: Consider the polynomial 


T(z) DH + Ed2), (2.16) 
where 
To(z) = 2" + B "+--+ Bz On Ken en (2.17) 


and £,(z) is a polynomial of any degree such that lim E,(z) = 0 as € — 0. 
Let the roots of To(z) be denoted by 21, 22, „ zn- Then J. (z) has at least 
n roots 21 (e), , Zn(€) such that 


limz,(€) = ze, &k=1,2,---,n. (2.18) 


By the lemma, the non-zero roots of TÏ) (w) approach those of T? (w) 
as € — 0. The total number of non-zero roots of all the TË Vs is therefore 


(nı - 0) +(n2—m) T (n nm-1) =n. Q.E.D. 
In our example (2.7), all of the non-zero roots of the associated poly- 


nomials TË Xw) appear to be regular. Starting with (2.10), we note, by 
inspection, that its regular roots are of the form 


w = do + age® + O(e°). (2.19) 
Substituted into (2.10), (2.19) implies that 
1+a3=0 

Zagas + a5 = 0. (2.20) 


Thus 
40 = (21) 07 = {ei*/5, ei, 5/3} 


as = 43/3 = 40/3. (2.21) 
The first three roots of (2.7) are therefore of the form 
25 (€) = wa (€) = e191 + / + Ofe7)], E=1,2,3. (2.22) 
The regular, non-zero roots of (2.11) are, by inspection, of the form 
w = bo + bie + O(¢?). (2.23) 
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This expansion, substituted into (2.11), implies that 
b3 + b§ = 0 (2.24a) 
362d, + 6b5b, + 250 = 0. (2.24b) 
The non-zero roots of (2.24a) are, again, the cube roots of —1: 
bo = (—1)¥3 = {eit/3, er, e5713}, (2.25a) 
With 63 = —1, (2.24b) yields 
b = —(2/3)b5. (2.25b) 
Thus the next three roots of (2.7) are of the form 


21 (e) = e~? we (€) = gy = (2x3) ce /s. O(e?)], 
k=4,5,6. (2.26) 


Equation (2.12) was the result of the intersection of three lines on the 


graph (see Fig. 2.2). The number of terms in TY Xw) is equal to the number 
of intersecting lines. When more than two lines intersect, we must consider 
the possibility of TË Xw) having multiple roots. If multiple roots occur, 
as happens in (2.12), then the order of the correction terms is somewhat 
harder to deduce. 


Exercise 2.3. By finding exact solutions, deduce the form of 
the perturbation expansions for roots of 

(a) 22 ＋ 425 -e 

(b) 22 ＋722 11e 


What interesting differences do you see in the two cases. 


Exercise 2.4. Compute the dominant term and one cor- 
rection term for (2.12) and (2.14). You will have to figure the 
appropriate order of the correction term based on experience 
gained in Exercise 2.3. 
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In general, the non-zero roots of the polynomials T) defined by 
(2.14) need not be regular: the a’s in (2.3) and the associated proper values 
and exponents, (pj, ej), may be non-integer rationals or TY) (u) may have 
repeated roots. Thus to obtain regular expansion, new parameters must be 
introduced. 


Let 
6 = 605, (2.27) 
where q j is the least common denominator of the set of exponents (O,. , an- 
np; }. Then from (2.14) 
RY(w) = T4)(w; 8%) = ee P(B%?s w; BY), (2.28) 


where Tw) = TO; e) and P,(z) = P(z;€). The roots of TY) (w) are 
identical to those of uy (w), but the non-zero roots of the latter will have 


regular expansions in G of the form 
w(B) = bo +b +--+ bnp + O(N). (2.29) 
In summary: 


Every root of the polynomial (2.1) is of the form (2.2). The set 
of exponents 2.6 determines a set {p1, , Pm} of proper values. 
With each proper value one introduces a new parameter ĝ via 


(2.27) and an associated polynomial RY Xw) defined by (2.28). 


The simple non-zero roots of RY (w) have regular perturbation 
expansions in f. The total number of non-zero roots of all the 


RS (v) is n. These yield expansions for each of the roots of 
(2.1). 


To pull together all the ideas in the chapter, let us construct expansions 
for the roots of the polynomial 


P.(z) = 1 — € + €(2 + 362) 2 — (16 — €)z4 ＋ 664 6028. (2.30) 


In outline, proceed as follows: 


1. Set z = ew in (2.30) and determine the set of exponents: 


Pe(ePw) = 1— e + ett (2 + 36% - e (16 — e) 
+e7*9P(4 - e + e) (2.31) 


E={0, 1+p, 3 + 4p, 2+6p}. (2.32) 
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2. Determine the pairs (pj, ej) of proper values and minimal exponents 


(with the aid of graphs or otherwise) and list the associated polyno- 
mials TY) (w) : 
(3/4, 0) : TO) (w) = 1 — e + e7/4(2 Ze) — (16 - e). 
7 / (4 ¢ +) w® (2.33) 
(5/2, —13) : T)(w) = €19(1 — €) + e??? (2 + 3?)w 
—(16—e)wt+(4—e+e)w® (2.34) 


For each j, determine gj, set € = 6%, and list the associated polyno- 


mial RY(w): 
e=: Ro (w) =1-p4+A"(2+ 36° )w 
—(16 — 6*)w* + 67°(4 — p4 + i (2.35) 
=f: RY(w) =P- G) + pP + 364)w 
—(16 — H) w + (4-6? + p'u. (2.36) 


Each R (v) has non-zero roots of the form (2.29). Substitute this 


expression into the identity RÖ )(w(B)) = 0, collect and equate to 
zero coefficients of like powers of 5, and solve, one-by-one, for the 
unknown coefficients bo, b1, ---- 


From (2.36), the roots of Rv) will have the form 

w(B) = bo + bapt + O(67). (2.37) 
However, from (2.36), the roots of NG (v) will have the form 

w(8) = bo + b28? + 0(6*). (2.38) 


The reason that (2.38) does not have an O(*) error term is that 
1006) must be raised to the fourth power which will produce terms of 
the order shown in (2.38). 


From (2.36) 
1— 1664 =0 


— 1 — 64b3b4 + b4 = 0. (2.39) 


Hence 1 
bo = (1/16)1/4 = * * 1, 2, 3,4 
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b4 = —(15/64)bo. (2.40) 
And from (2.36), 
—16b4 + 4$ = 0 
— 64b3b2 + 56 + 245662 — of = 0. (2.41) 
Hence 
bo = (4)'/2 = 2(-1)', k=5, 6 
b2 = (3/32)bo. (2.42) 
5. Write down the roots 21 (6), , ze(€) of (2.32) via the change of 


variable z = Pw: 
ze(€) = 1/2e-**/21] — (15/64)e+ O(e7/*), k= 1, 2, 3, 4 (2.43) 
ze (e) = 2(—1)*e7*/7[1 + (3/32)e + O(e?)], K = 5, 6. (2.44) 


Exercise 2.5. Compute explicitly the O(e2) term in (2.44). 


Exercise 2.6. Compute the first two non-vanishing terms in 
the expansions of each of the roots of 


(a) IT + etz? ＋ 25 
(b) 1— 22 722 ＋T E25 
(c) 1722 +2? + e424 ＋ 6725. 


Exercise 2.7. The governing equations of the Morley- Koiter 
theory of circular cylindrical shells can be reduced to a single 
eighth order partial differential equation. Separation of vari- 
ables yields a constant coefficient ODE whose associated poly- 
nomial is 
* = m?)? (z? m $ 1)? $ 2 

where m = 0, 1, 2, --- and é is proportional to the thickness 
to radius ratio of the shell. Determine one term expansions with 
remainder estimates (i.e. O-terms) for each of the roots, taking 
note of the special cases m = 0,1. [Hint: The polynomial is of 
the form a? +b? = (a +ib)(a—ib). As the roots of a polynomial 
with real coefficients must occur in complex-conjugate pairs, the 
roots of a- ib will be the conjugates of the roots of a+ ib. Also, 
each factor is quadratic in 2°]. 


Chapter 3 


Singular Perturbations in 
Ordinary Differential 
Equations 


The simplest type of ODE is linear and has constant coefficients: finding its 
general solution hinges on finding the roots of the associated polynomial. 
If the coefficients in the DE depend on a parameter e then so do the roots 
of this polynomial. In general, a study of the behavior of these roots is not 
sufficient to infer the behavior of the solution of the DE because 


1. A non-homogeneous term in the DE can give rise to a term in the 
general solution whose behavior depends only, in part, on the roots 
of the associated polynomial. 


2. The imposition of initial conditions (IC’s) or boundary conditions 
(BC’s) may result in the appearance of e in the constants multiplying 
the homogeneous terms in the general solution. 


3. The domain of interest may depend on e. 


4. Most importantly, the solution of the DE will be a function of both 
e and the independent variable, say t, which immediately raises the 
question of whether certain approximate solution are uniformly accu- 
rate for all values of t in the domain of interest. 


Fortunately, when € is small, characterizing the e-dependence of the 
solution of a DE with IC’s or BC’s simplifies. The rest of this book is 
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devoted to this task. This chapter is intended to acquaint you with two 
simple problems that are typical of those we shall study later. 
First, we consider the initial value problem (IVP) 


Ley) = §+2y+y=0, lel<1, O ct, y= dy/dt (3.1) 


y(0) =0, (0) =1. 


If € is small, we might expect the solution of (3.1) to be approximated 
accurately by the solution of the reduced problem 


Lo(Y)=Y+Y=0, 0<t (3.2) 


Y(0)=0, Y(0)=1, 


which is 
Y(t) =sint. (3.3) 


Substituting (3.3) into the DE in (3.1), we see that, on the average (here 
meaning over one period, t = 27), the magnitude term 26 is uniformly 
small, of order e, compared to the magnitudes of each of the terms Y and 
Y. However, 


Small disturbances acting for long times can have large effects 
as we see upon comparing (3.3) with the exact solution of (3.1), 


e sin( V1 — ert) 


12 (3.4) 


y(t, e) = 
(See Fig. 3.1) 
Exercise 3.1. Verify (3.4). 


Because of the factor e~“ in (3.4), the amplitude of the exact solution 
either grows without bound (e < 0) or decays to zero (e > 0) as t — 
oo. Thus on the semi-infinite domain 0 < t, Y(t) is not an accurate 
approximation to y(t, €) and we say that problem (3.1) exhibits a singularity 
in the domain. Such behavior has no counterpart in the algebraic problems 
we studied in Chapter II, as the polynomial z? + 2ez + 1 associated with 
(3.1) is regular. 
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Figure 3.1: Solution of Eq. (3.1), y(t, €) (solid line), and Eq. (3.2), Y(t). 


On any finite domain, 0 < t < T, we may apply the mean value theorem 
to (3.4). With ye = dy/de, X. VI e and e, is some number between 
0 and e, we have 


y(t, €) = y(t, 0) + ye(t, e.) e 
t sin (A. t) e. t cos (A. t) 
» 902 
c. sin (A. t) 
* 
sint T Ole), O St ST. 


= sint + e~*** |- 


(3.5) 


Thus we conclude that (3.1) is a regular perturbation problemif0 < t < T, 
but not if 0 < t. How to tame singularities arising from (semi-)infinite 
domains will be studied in Chapter V. 
A second interesting phenomenon occurs in the boundary value problem 
(BVP) 
L.) Sch = , 0<t<1, (3.6) 


y(0) 0, 1) = 1. 
Immediately we see that we have a singular perturbation problem because 


the reduced equation A 
Lo(Y)=Y+Y=0 (3.7) 


is of first order. Its solution 


Y(t) Ae, A constant, (3.8) 
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cannot possibly meet the two boundary conditions in (3.6). Two questions 
now arise: Does (3.8) in any way approximate the exact solution of (3.6), 
and, if so, how can A be determined? 


Exercise 3.2. Find the exact solution of (3.6) 


We could answer our two questions by examining the solution to the 
preceding exercise. Instead, let us try to obtain, heuristically, an accurate 
approximate solution to see what it tells us. We do so to find an approach 
to more complicated problems where exact, closed-form solutions may not 
be available. 

The solution of the DE in (3.6) is the sum of two independent homo- 
geneous solutions. There is some hope that one of these might be closely 
approximated by (3.8) since (1), if we substitute (3.8) into the unreduced 
DE, the neglected term eV is uniformly small compared Y and Y ( i.e., 
compared to each of these terms individually, not compared to their sum, 
which is zero); and (2), the domain of the independent variable is finite. 
(We have deliberately considered a boundary value problem on a finite do- 
main to avoid compounding the problems of a singular domain with that 
of a singular equation. In some problems, both difficulties are present si- 
multaneously). 

Any second, independent solution of the DE in (3.6) must have the 
property that the first term ej, for some values of t, is of the same order 
of magnitude as one or more of the remaining terms; otherwise we would 
obtain again the reduced DE as a first approximation. 

The polynomial associated with the unreduced DE, 


P. (z) = ez? +241, (3.9) 
is singular. To study its singular root, we make the substitution z = e~!w 
and obtain 

eP,(e“!w) =w +w +e. (3.10) 


Making the equivalent substitution t = er in (3.6), we have, with y’ = 
dy/dr, 


eL(y) =y" +y +ey=0,0<7< E, y(0)=0, ye") = 1. (3.11) 


If we now seek a solution of the transformed problem (3.11) such that the 
term e is small compared to y” and /, we are led to consider the simplified 
DE 

w"+w'=0 (3.12) 
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Figure 3.2: Graphs of ee~* (solid line) and e(e — e /) with e = 0.05. 


whose general solution is 
BTCC = BIT Ce-. (3.13) 


The constant B must be set to zero, otherwise e will not be small com- 
pared to B“ and B’ (which are zero), as was assumed. 
Adding the two independent approximate solutions (3.8) and (3.13), we 
obtain 
y= Ae T Ce /, (3.14) 
which we hope will lead to a good approximation to the solution of (3.6). 
Imposition of the BC’s in (3.6) yields 


2 = = 
1 . te * 


= e(e~* e /) + O(e /). (3.15) 


In Fig. 3.2, we compare the graph of the first term in the second line of 
(3.15) with the graph of ee. 


Exercise 3.3. Use the result of Exercise 3.2 to show that 


vlt, €) ele e = Ole). 
(Hint: Try to arrange the solution of Exercise 3.2 so that it 
looks as much like (3.15) as possible. Then drop transcenden- 
tally small terms. These are terms that go to zero exponentially 
at every point in the interval.) 
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The term e / is called a boundary layer (or the boundary layer con- 
tribution to y(t, e)) because it is significant only in a narrow layer of width 
O(e) near the boundary t = 0. As setting € = 0 in (3.6) often corresponds to 
replacing a physical model by a simpler one ( e.g., a shell by a membrane), 
we shall say that such a BVP exhibits a singularity in the model. We can 
now answer, to some extent, the two questions we posed earlier. First, Ae! 
does approximate the solution of (3.6) outside a boundary layer at t = 0. 
For this reason, Ae~* is called the interior contribution to y(t, e). Second, 
because there is no boundary layer neart = 1, A can be determined from 
the BC att = 1. 

In the chapters to follow, we shall develop systematic methods for re- 
ducing singular perturbation problems to regular ones. These will then be 
solved by constructing regular perturbation expansions. 


Exercise 3.4. Explain how one or more of the points dis- 
cussed at the beginning of the Chapter is illustrated by each 
of the following problems. Also explain why each is a singular 
perturbation problem. (Note that all may be solved exactly.) 


(a) y+ (1+ ey=0, Oct = 
y(0)=0, ym) =e, p>. 

(b) ty ＋ ty = O, e<t<l 
ye) = 1, y(1) =0. 

(c) y+ty=ecost, 0<t 
y(0)=0, 500) S 1. 


Chapter 4 


Periodic Solutions of the 
Simplest Nonlinear 
Differential Equations. 
Poincaré’s Method 


A Physical Model. Consider a frictionless cart of mass m attached to a 
spring that is anchored to a rigid wall (Fig. 4.1a). Let c denote the un- 
stretched length of the spring and let z denote an additional displacement. 
From experiments on the spring we can construct a force-displacement 
curve. If the curve has a non-zero slope E at zero displacement, then 
it can be presented as a dimensionless stress-strain curve, as indicated in 
Fig. 4.1b. 


Assume that the stress-strain relation is valid for any two springs having 
the same material properties and shape and let values of the stress F/E 
at a strain z/c be denoted by f(z/c). Further, assume that tension in the 
spring produces extension and that compression produces contraction, i. e., 
assume that (z/c)f(z/c) > 0, z # 0. The equation of motion of the cart 
then takes the form 

dz 


my = -F = -Ef(z/c). (4.1) 
To be specific we take as initial conditions 
—— ei, (4.2) 


dT 
45 
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Figure 4.1: (a) Mass attached to Figure 4.1: (b) Stress-strain rela- 
a nonlinear spring. tion for the spring. 


We wish to study the motion of the cart, especially when the initial 
strain, zo/c, is small. To this end it is convenient to nondimensionalize our 
IVP. A natural choice for a dimensionless displacement is æ /o; another is 
the strain itself, z/c. The first choice makes the IC’s parameter free; the 
second, the DE. Because f(z/c) is essentially arbitrary at this stage, we 
shall take the strain as the dimensionless measure of displacement. Alto- 
gether then, if we set 


z=cu, 20 S ce, TS t Vm / E, (4.3) 
our IVP takes the form 
ü+ F(u) =0, O St, u(0) e, u(0) = 0, (4.4) 


where ù = du/ dt. Note that the IVP (1.87) for the large amplitude motion 
of a pendulum corresponds to taking f(u) = sin u. 

Before attempting to construct a suitable representation for the solution 
of (4.4), let us consider the simpler problem of determining, as a function 
of e, the (dimensionless ) period of oscillation of the motion. First, though, 
we should assure ourselves that (4.4) has a periodic solution. 

We begin by noting that if we multiply both sides of the DE in (4.4) by 
ù, the resulting equation can be written 


K+V=0, (4.5) 


where 
=-u (4.6) 


denotes the (dimensionless) kinetic energy of the cart and 


V.. I "njü (4.7) 
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denotes the (dimensionless) potential energy of the spring. Thus (4.5) im- 
plies the equation of conservation of energy, 


K +V =C, a constant. (4.8) 


From the IC’s in (4.4), C = V (e€). Solving (4.8) for ù, we have 


ú =+/21V( — Vu). (4.9) 


To exploit our physical intuition, it is useful to recognize that (4.9) 
also represents the equation of motion of a bead that starts at rest from a 
height V (e) and slides without friction along a wire of height V (u). Here u 
measures distance along the wire, as indicated in Fig. 4.2. 


Figure 4.2: Bead Sliding along a wire whose height is a function of distance 
along the wire. 


As Fig. 4.2 suggests, a necessary condition for periodic motion is that 
there exist a constant —e, such that V(—e.) = V (e). For if the elevation of 
the wire to the left of u = 0 never reached V(e), the bead, once in motion, 
would continue to move forever to the left. (If the constant —e, exists, it 
must be unique. This follows from (4.7) and the fact that uf (u) > O, u # 0, 
which together imply that V is a strictly monotonically increasing function 
of |u]; i.e., V(u) is 1:1 on u > 0 and u < 0.) 

To extract further information from (4.9), we introduce the very useful 
concept of a phase plane, the set of all ordered pairs of real numbers (u, ù). 
The phase portrait of the solution of (4.4) is the graph of the two functions 
defined by taking the + or — sign in (4.9). Since V(e) — V(u) steadily 
decreases from V (e) to 0 as u increases {decreases} from 0 to e{—e.}, the 
graphs of these two functions form a single closed curve, as depicted in Fig. 
4.3. 

Equation (4.9) does not tell us how a point on the phase portrait varies 
with t. For this information, we return to the DE in (4.4) and write it in 
the form 

(u) = —f(u). (4.10) 
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Figure 4.3: Phase portrait associated with K + V = V(e). 


This equation and Fig. 4.1b imply that if u is positive, then ù is a decreasing 
function of time, and vice-versa. Thus a point on the phase portrait moves 
clockwise, as indicated in Fig. 4.3. 

The period of oscillation P is the time of one traverse. From (4.9) and 
Fig. 4.3, : P 

u 
p Oa 

At the limits of integration the integrand has singularities. These must be 
integrable if the motion is to be periodic. 

(To those who know some complex variable theory, we point out that if 
we regard u as complex, cut the u-plane from —e, to €, assume that V (u) 
can be extended analytically to a neighborhood of this cut, and choose the 
branch of the square root in (4.9) that reduces to ,/2V(e) when u = 0, we 
obtain the elegant formula 


P= f — h: 
c V2[V(e)- V(u)] 


where the integral is taken in a clockwise sense along a loop C surrounding 
the cut from - e. to €, as shown in Fig. 4.4.) 

Computation of P(e) for |e] < 1 Directly by Formula. To be 
specific, we assume that 


(4.11) 


(4.12) 


f(u) Su 55 u? <6.) (4.13) 


17 (u) is the two term Taylor polynomial for ain u about u = 0. 
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Figure 4.4: Contour for the complex line integral in Eq. (4.12). 


From (4.7) the potential energy is 


1 u? ut 
V(u) = f 6-30 4 (4.14) 
In this case €, = e, and (4.11) reduces to 
du 
PC we uU Hef, ( 
To facilitate integration let 
u= esing, O < a < 7/2. (4.16) 


Then 


1... 
7 f {1 —(e?/12)(1 + sin? a)}4/2 
17/2 
= af [1 + (c?/24)(1+sin?a) T. Ida, e C 6 
0 
= In[1 + 2/16 ＋ O( e)], e geg <6. (4.17) 


Thus, the larger the initial amplitude é, the longer the period of oscillation. 
This is what we expect physically, since (4.13) represents a “soft” spring, 
i.e. a spring whose stiffness decreases as it stretches. 

Computation of Pe) for |d & 1 Indirectly by Poincaré’s Method 
The solution (4.17) depended upon our being able to partially integrate the 
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DE in (4.4). Now we consider an alternatative method that does not re- 
quire this. In addition to obtaining an approximate formula for P(e), we 
shall obtain a good approximation to the motion itself. 

With (4.13) and the change of variable 


u=ev, P= 8, (4.18) 
(4.4) reduces to 


ü+v- 9005 =0, 0 St, v(0) =1, 500) =0. (4.19) 


For a wide class of IVP's that includes (4.19), there is a theorem? that 
says the solution is analytic in G for t and 2 sufficiently small. Thus we are 
guaranteed the regular expansion 


v(t, B) = volt) + pvi (t) + +--+ BY w(t) + NRG, A), (4.20) 
where, for some T > 0, 
Rysilt, 6) ON), KT. (4.21) 


This useful result reduces the search for a function v(t, H) of tuo variables 
to a sequence of searches for functions vo(t), vi (t). .. of one variable. We 
will not use the Theorem because its proof requires too many new tools, and 
because it may fail to apply in more complicated problems where physical 
experience suggests that a representation of the form (4.20) is reasonable. 
Instead, we shall simply assume that (4.20) and (4.21) hold and work out 
the consequences. Those who are unhappy with this attitude may regard 
all that follows as merely a formal way of obtaining a conjectured form of 
solution which then is to be justified rigorously. 


Why the Unmodified Regular Expansion Won’t Work. To pro- 
ceed, we must assume that the derivatives ù and ö also have expansions 
with remainder estimates of the form (4.20) and (4.21). Substituting these 
expansions along with (4.20) into (4.19), we have 

(ïo + Bin +-+ BY in + Ën) + (v0 + pvi ++ BY ow 

NR) - (1/6) (vo + Bu + +--+ BX uy Rr) S O, (4.22a) 
vo(0) + 8v1 (0) + --- + BY uy (0) + Rv41(0, G = 1, (4.22b) 
ùo(0) + Bùi (0) +--+ + BY ùn (0) RN, 8) = 0. (4.22c) 


2 Coddington & Levinson, Theory of Ordinary Differential Equations, McGraw-Hill, 
1955, Theorem 8.4, pp. 36-37. 
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In (4.22a) we expand the cube in the second line and then collect coefficients 
of like powers of @ to obtain 
üo+vo + Ali, + vı — (1/6)v3] 
+ B? [jig + v2 - (1/2) vil + ---+O(6%t4)=0. (4.23) 


By the fundamental theorem of perturbation theory, the coefficient of 
each power of p in (4.23) and the IC’s of (4.22b) and (4.22c) must vanish. 
This yields the sequence of IVP’s 


Üo + vo = 0, vo(0) =1, vo(0) = 0 (4.24a) 
51 + v1 = (1/6)vp , (0) = (0) =0 (4.24b) 
üz + v2 = (1/2)vĝvı, v2(0) = ù2(0) = 0, (4.24c) 
etc. 
The solution of (4.24a) is 
vo (t) = cost. (4.25) 


Before proceeding farther, let us see how well vo(t) approximates v(t, H). If 
the expansion (4.20) and remainder estimate (4.21) are valid, then for any 


fized T > 0, 
u(t, 8) — vo(t)| = O(8), ltl < T. (4.26) 


The period of vo(t) is 2m whereas (4.17) shows that the period of v(t, G) 
depends on B. Thus, as t increases, vo(t) and v(t, G) grow slowly out of 
phase (Fig. 4.5). 

Indeed, ift = O(@-") then, by (4.17), all we can assert is that |v(t, 2) — 
vo(t)| = O(1). Clearly we are dealing with a perturbation problem in which 
there is a singularity in the domain. Will taking a two term approximation 
to v(t, G) improve things? 

Inserting (4.25) into the DE in (4.24b) , we obtain 

cost cos at 


r (4.27) 


upon using the trigonometric identity 
cos*t = Š cost + F cost. (4.28) 


The solution of (4.27) that meets the IC’s in (4.24b) is easily found to be 


tsint cost —cos*t 
v(t) = 16 a (4.29) 
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Figure 4.5: Graphs of v(t, £) (solid line) and vo(t) showing how the two 
functions drift out of phase. 


For any fixed T > 0, (4.21) tells us that 
lolt, A) — v(t) — Hunde Ode. (4.30) 


However, it is clear from (4.29) that on the full interval of interest, 0 < 
t, volt) + v(t) is a poorer approximation to v(t, B) than vo(t) alone. The 
culprit is t sin t which, unlike the exact solution, is non-periodic and grows 
without bound as t — oo. (Such a term is called a secular term.) More- 
over, the approximation vo (t) + fv; (t) fails to reflect in any simple way the 
variation with 2 of the period of oscillation. 


Exercise 4.1. Compute the period P by observing that P 
equals 4 times the value of t at which v(t) + v: (t) first crosses 
the t-axis. Compare results with (4.17). 


If we interpret (4.27) as the equation for the forced motion of a linear, 
undamped oscillator, then the factor t sint in (4.29) represents the response 
of the oscillator to a resonant driving force, (1/8) cos t. 

In summary 


Though the solution of the IVP (4.19) admits a regular expan- 
sion, the remainder Ri (G, t) is not uniformly O(ØN+!) for 
0 < t. To get an approximation to the exact solution, useful 
over times long compared to the period of oscillation , we must 
seek another type of expansion. 
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Poincaré’s Method of Strained Coordinates. Let us reconsider 
our first approximation to v(t, f). The function vo(t) = cost agrees with 
v(t, H) in amplitude, but has a slightly shorter period. It is this discrepancy 
that drives vo(t) and v(t, £) slowly apart. To resolve the difference suppose 
that we replace the first approximation by cos At, where À is a constant. 
Then there must exist a value of A such that cos At and v(t, G) have iden- 
tical periods. Poincaré’s method is a systematic procedure for determining 
successively more accurate approximations to À. 

The variable 

T= (4.31) 
is called a strained variable. The first step in Poincaré’s method is to 
introduce T as the new independent variable. Then (4.19) takes the form 


3 =p” =0, 0< T, 00) 1, % %% (4.32) 
where v' = dv/dT. Next, because À — 1 as 6 — 0, it is natural to assume 
that A has a regular expansion of the form 

A(B) = 1+ NM + B?A2 + +--+ BN AN + ONY), (4.33) 


where the unknown constants Al, Az, . are to be found. Finally, we assume 
that the solution of (4.42) itself can be represented in the form 


(T, B) = zo(T) + Bai (T) + G +---+ BX zw RN, G), (4.34) 


where 

Ri = O(6%*"), (4.35) 
uniformly in T. Substituting (4.33) and (4.34) into (4.32) and equating 
to zero the coefficients of successive powers of H, we obtain the following 
sequence of DE’s and IC’s. 


60: 20 ＋ 20 = 0, 2(0)=1 24(0)=0 (4.36a) 
B:ata= ae — 23120, 21(0) = z1 (0) = 0, (4.36b) 


BP: 4 + za 821 f — (AP + Da), 
z2(0) = 24(0) =0, (4.36c) 


etc. 
The solution of (4.36a) is 


zo(T) = cosT. (4.37) 
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Substituting this into (4.36b) and using the trigonometric identity (4.28), 
we obtain 

21 + 21 = (1/8 + 2A1) cos T + (1/24) cos 3T. (4.38) 
The cos T-term on the right hand side of (4.38) —a resonance-producing 
term—will give rise to a term proportional to T'sin T in the solution for 
21. Recall that it was such a term in our earlier solution for vi (t) that 
made vo(t) + Hui (t) such a poor approximation to v(t, B). However, now 
we have the constant Ai at our disposal, and we shall choose it to suppress 
the resonance-producing term in (4.38). That is, we take 


1 


* ig (4.39) 
With Ai in hand, our first approximation to v(T, H) takes the form 
zo(T) = cos T = cos At = cos[1 — (1/16)8 + 007) t. (4.40) 


The period P of oscillation of 20 [which we have required be equal to that 
of v(t, B)] satisfies the relation 


[1 - (1/16) + 005% P = 2r, (4.41) 
Be, 

P = 2r[1 + (1/16) + 0(8")], (4.42) 
which agrees exactly with (4.17). 


Now comes an extremely important observation. If we wish to use the 
results we have obtained so far, namely zo(T) and Al, and no more, the 
best approximation we could use for v(t, G would be 


v(t, 8) & cos[1 — (1/16) ft. (4.43) 


Recall that the first approximation solution, cost, gave an error of O(£) for 
t = O(1), but a larger error of O(1) for t = O(@~"). However, (4.43) now 
gives an error of only O(£) fort = O(@-1) and does not lead to an error 
of O(1) until t = O(@-?). 

It seems clear what to expect as we compute higher and higher approx- 
imations to v(t, G). The right hand sides of each of the DE’s in (4.36) will 
contain a resonance-producing term which must be suppressed by a proper 
choice of A2, A3, . In this way we obtain successively better approxima- 
tions to the period P. At the same time, the knowledge of z2(T), z3(T), --- 
allows the motion itself to be computed more and more accurately. 

For example, if we compute 21 (T) but not Az, we can expect zo(T) + 
521 (T) to approximate v(t, H) to within an error of O( g 2) for t = O(@~"). 
If we also determine Az, we obtain the same degree of approximation to 
v(t, G) but for t = O(8-?), etc. 
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Exercise 4.2. Compute zı(T) and A2. Use your value of Az 
to compute the next term in the expansion for the period given 
by (4.17). Check your result by direct computation from the 
second line of (4.17). 


In Appendix B we prove that the remainder in (4.34) is uniformly small 
for all T > 0. That is, (4.34) is not an assumption but a consequence of 
(4.31) to (4.33). 

In summary 

By partially accounting for the effect of the small parameter 
hin the first approximation to v(t, G), i.e., by replacing t by 
T = \(@)t, we have reduced a singular perturbation problem to 
a regular one. However, there is a limitation. If only Ai, , Ax 
have been computed, then (4.34) can be used to approximate 
v(t, H) only if t = O(@-*). Since |8| < 1, this may be no real 
limitation in practice. 


Exercise 4.3. Consider the IVP 
ü TuT cu? = 0, 0<t, u(0)=1, u(0) =0. 

(a) For what range of values of e do you expect periodic 
solutions to exist? 

(b) If «| < 1, try to solve for u(t, e) by assuming a regular 
perturbation expansions in powers of e. Carry your 
calculations far enough to indicate where and why the 
assumed expansion fails. 

(c) Apply Poincaré’s method to avoid the difficulties you 
found in (b). 


Exercise 4.4. Use Poincaré’s method to construct the first 
term in the regular expansion of the solution to the IVP 


1 1 à 
ü+v+t 25ßv =0, KI, 0 St, v(0) = 1, 50) = 0. 


Hint: cos t l cos t is an even function of period 27 and hence can 
be represented by the Fourier series Co a, cos(nt), where 


ao = (1/7) 7 cos t|dt 


an = (2/7) [ cos(nt) cost |cost|dt, n > 1. 
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Forced Motion of a Nonlinear Oscillator. Suppose that the fric- 
tionless spring-mass system of Fig. 4.la is subject to an external, time- 
dependent force. Its dimensionless equation of motion can then be written 


t+ f(u) =g(t), 0<t. (4.44) 


If f(u) = u, i.e., if the DE is linear, then the solution of (4.44) is the 
sum of an homogeneous solution plus a particular solution: 


u=un(t) + w(t) 
= u(0) cost + ù(0) sin t + Z sin(t — r)g(r)dr. (4.45) 


The last term on the right of (4.45) is called the forced response. 


Exercise 4.5. Verify that (4.45) satisfies (4.44) if f(u) = u and 
that up(0) = u,(0) = 0. 


If, for a general f, we set f(u) = u+([f(u) — u] in (4.44) and take the 
term in brackets to the right side, then (4.45) becomes a nonlinear integral 
equation for u: 


u = u(0) cost + u(0) sin t 
+ | sin(t — r)[u(r) — f(u(r)) + g(r) ]dr = (u). (4.46) 


Assuming |u| < 1, we might attempt to obtain a sequence of approximate 
solutions to (4.46) by the iteration process un gi = K(un), uo = 0. Alterna- 
tively, we might ask: can Poincaré’s method be used to obtain approximate 
solutions? In general the answer is no, but in an important special case, 
that we shall now discuss, the answer is yes. 

If the dimensionless, external force is harmonic, say 


g(t) = esinut, (4.47) 
then if f(u) = u, (4.45) yields 


wsint —sinwt 
2 — 1 


— 2 ban: t cost) asw 1. (4.48) 


up = 
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Thus, if w Æ 1, the forced response of the linear spring-mass system is har- 
monic, but ifw = 1, the forced response is non-harmonic with an amplitude 
that grows without bound as t — oo. Such resonance is of great practical 
concern in elastic structures. 

Invariably, if an elastic structure undergoes large, non-rigid deforma- 
tions, nonlinear internal forces appear. Thus, in the simplest of all elastic 
models, the spring-mass system, it is natural to ask what happens to solu- 
tions of 

ü + f(u) = esinwt (4.49) 
ifw 1. To be concrete, let us assume that f(u) is odd and has a repre- 
sentation for small u of the form? 


f(u) = u+ ku T Ou), k 00). (4.50) 


If w = 1 and if at some instant |u| < 1 we expect the resonant term 
et cost to begin to dominate the motion. But as |u| grows, the nonlinear 
term kus in (4.50) will come into play, moving the natural frequency of the 
force-free system away from 1. Thus ful t; v, e)] need not approach oo as 
wl. 

For the nonlinear term ku? to counteract the resonance-producing effect 
of the forcing term esinwt, we must have u = O(e!/5). This suggests the 
change of variable 


u = (/ E) J, 6 = ke.. (4.51) 


Futthermore, as we wish to determine u not just at w = 1 but ina neigh- 
borhood of w = 1, we set 


w =1+4 fo, Tut. (4.52) 


The parameter o may be called the detuning. Substituting (4.50) to (4.52) 
into (4.49) and dividing by (8/|k|)!/?, we obtain“ 


(1+ Bo) +v + sgn(k)bv? = BsinT + 0(8*), v' =dv/dT. (4.53) 


To obtain approximate solutions to (4.53) for arbitrary IC’s consistent 
with the assumption that v = O(1), we need the two-scale method discussed 
in the next chapter. However, if we are content with any solution of (4.53) 
that is of period 2—take what IC’s many come then Poincaré’s method 
suffices. 

If f(u) = u kus, (4.49) is called Duffing’s equation. 
4sgn(z) is the signum or sign function: sgn(z) = 1,2 > d sgn(0) = 0, sgn(x) = 
—1, 1 <0. 
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Thus, we look for solutions of (4.53), uniformly valid in T, of the form 
oT; 8) = vo(T) + Bu (T) + O(6"), (4.54) 


the dependence on e being understood. Substituting (4.54) into (4.53) and 
equating to zero coefficients of 8° and H we obtain 


to + =0 (4.55a) 
v” +v = sin T- ovo” — sgn(k)vg. (4.55b) 


Substituting the solution of (4.55a), 
vo = Ao cos T + Bo sin T, (4.56) 


into (4.55b) and using the trigonometric identity (4.28) plus 


cos? t sin t = isint + 4 sin 3t (4.57) 
— 1 
cos t sinꝰ t = zt- 4 cos 3 (4.58) 
3 (4.59) 
4 4 3 : 


we obtain 
vf +01 = Aolo — 4048 + Nee 
+[1 + 9 ⁰ = Bo( AB + Bĝ)]sin T 
+7 Aol 383 — A?) cos 37 — 4543 — B2) sin 3T. (4.60) 


To avoid resonant terms in vi, the coefficients of cos T and sin T must be 
set to zero. This yields 


40 = 0, 1+ 00 — 455 =p; (4.61) 
Fig. 4.6 is a graph of Bo vs. o (computed by regarding ø as a function of 


Bo). Note that, depending on the value of ø compared with a, = (3/2)*/3 = 
1.717 - - there may be one, two, or three solutions to (4.49) of the form 


u = (8/|k|)!/?[Bo(c) sin wt + O()]. (4.62) 
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Figure 4.6: Amplitude of the harmonic solutions of Duffing’s equation as a 
function of the detuning. 


We have taken only the first step in the analysis of Duffing’s equation. 
The treatment of arbitrary IC’s, the inclusion of damping, and the elucida- 
tion of the behavior of the oscillator in a neighborhood of ø = o, may be 
found in more detailed works on perturbation theory.® 


Exercise 4.6. Suppose that w 3 and that k < 0. Set w = 
3(1+ 80), T =(1+c)t and show that there exist solutions of 
(4.49) of the form 
E 1720 8in( 
u = (6/\k))"7[2sin( ut) + O(P) 


These are called sub-harmonic solutions. 


5See the two books by Nayfeh, Introduction to Perturbation Methods, Wiley, 1981 and 
Perturbation Methods, Wiley, 1973. 


Chapter 5 


Introduction to the 
Two- Scale Method 


The Damped Linear Oscillator. Assume that the relative velocity of 
a pendulum swinging in air is small enough for the drag on the bob to 
be viscous and proportional to the velocity. If the only other damping is 
assumed to come from a torque at the pivot proportional to the angular 
velocity, then the equation of motion is 


do j dð 
ma = = mg sin — ef. (5.1) 


where m is the mass of the bob, l is the length of the pendulum, @ is its 
angular displacement from the vertical, g is the gravitational constant, c is 
a damping factor, and t, is the time. Suppose that the pendulum, hanging 
at rest, is struck by a hammer at t, = 0. Then the IC’s are 


6=0, d0/dt.=w, (5.2) 


where mlw is the impulse imparted by the hammer. 
To work with the simplest model possible, assume that [t.) ] < 1, 
replace sin by d and nondimensionalize as follows: 


= tV1/g, 2e= c/(my/ig), 9 v0 s. (5.3) 
Then, with () = d )/dt, the DE (5.1) and IC’s (5.2) reduce to the IVP 
t+2v+v=0, 0<t, v0) =0, 0(0)=1. (5.4) 


61 


62 Introduction to the Two-Scale Method 


Exercise 5.1. Show that the exact solution of (5.4) is 


sin(t VI — e2). (5.5) 


e et 


a 


Why an Unmodified Regular Expansion is No Good. Ignoring 
for the moment Exercise 5.1, let us see what follows if we assume that 1) 
the solution to (5.4) has the regular expansion 


vlt, e) = 


v(t, e) = volt) + evi(t) T.. + eV un(t) + RAI (t, e), (5.6) 
that 2) tio, bi, , iv, and Ry41 exist, and that 3) 
Ry, RAI, RI = O(t). (5.7) 


Substituting (5.6) into (5.4) and invoking The Fundamental Theorem of 
Perturbation Theory , we get the following sequence of IVP’s: 


©: 50 ＋ vo =0, vo(0) = O, 5000) =1 (5. 8a) 
el: ü +v ＋ 2 S 0, (0) = 100) =0, (5. 8b) 
etc. The solution of (5. 8a) is 
volt) = sin t. (5.9) 
Hence (5. 8b) reduces to 
üi +v =—2cost, vi (0) = 10) =0, (5.10) 
whose solution is easily found to be 
v(t) = —tsint. (5.11) 


So far, we have 

v(t, e) = sin t - et sin t + R,(t, €). (5.12) 
Taking a peek at Exercise 5.1, we see that (5.12) is simply a Taylor expan- 
sion in € of the exact solution. 


Exercise 5.2. For t fixed and |e] < co < 1, show that 
|Ra(t,€)| < Ke. Note, however, that the term t sin t grows 
without bound as t — oo. Hence Rz is not uniformly small in 
t. 
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Motivating the Two-Scale Method. To obtain a regular expansion 
that is uniformly valid we must, somehow, pull a part of the effect of € into 
the first approximation to v. This is what the two-scale method of Cole 
and Kevorkian does. To provide physical motivation, let us introduce the 
dimensionless energy 


E= 20 +02). (5.13) 
Then with the aid of the DE in (5.4) we have 
E = Ae. (5.14) 


The IC’s in (5.4) yield ECO) = 1/2, so, upon integrating both sides of (5.14), 
we get 


E(t) = 5 —2 f "Pas. (5.15) 


Hence, E decreases monotonically. 

Let us estimate the relative decay of E over one period. To within an 
error of order e, we may replace ù by ùo = cost on the right side of (5.15). 
Thus 


an a Ae | — (5.16) 
0 
or, roughly, 
as x —2eAt. (5.17) 


Integrating and setting E(0) = 1/2, we get 
E(t) & 20 (5.18) 


Thus E decays to roughly 1/e of its initial value over a time t = O(e7*). 

This analysis shows that the solution for the damped linear oscillator 
contains two time scales, the period of oscillation, t * 24 = O(1), and the 
energy “half life“, t (26) 1 = O(e~1). These observations suggest that 
it may be fruitful to think of the exact solution v as depending on two 
independent variables, or two scales, t and 


r= et. (5.19) 


In the terminology of Cole and Kevorkian ,t is the “fast” variable and 7 
the “slow” variable. 
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From Ordinary to Partial Differential Equation. Assuming v = 
v(t, 7,€), we have by the chain rule and (5.19), 


_ dv _ dv, dvdr 
= "dt ot ardt 
dt or 
= vt + evr. (5.20) 
Likewise, 
d 
v= 45 0 = (v: + eur): + (v: + eur) 47 
= Vt + 2ever + 9 . (5.21) 


Inserting (5.20) and (5.21) into the original IVP, (5.4), we obtain 


ve + v + 2e(vir T vt) + eur + 207) = 0, 0<t, 7, (5.22) 
v(0,0,€) =0, vr (0, O, e) T evr (0, O, e) =1. 
(Things get worse before they get better!) 


The Modified Regular Expansion. Now assume that v has the 
uniformly valid expansion 


vlt, 7, e) 9 (t, 1) ＋ e ò (t. 7) +---+e% ` (t, 7) + Rvyilt,7, e), 
RN+ı = O(eN * (5.23) 


This expression, inserted into (5. 22), implies the following sequence of par- 
tial differential equations (PDE’s) and IC’s: 


eè : We + 0=0, 9 (0, 0) =0, v: (0,0) =1 (5.24a) 


el i Ui +0 +2(0er+ v) = o, v (0,0) =0, % (0,0)+ 8, (0,0) = 0, 
(5.24b) 


etc. 
Comparing the DE in (5.24b) to that in (5.8b) we see that the underlined 
term is new. It is this term which will give us the necessary freedom to 


suppress secular terms in ò. 


Exercise 5.3. Find the next equation in the (5.24) sequence 
by setting the coefficient of €? to zero. 


A First Look At Perturbation Theory 65 


Solutions. The PDE in (5.24a) looks, formally, like the ODE in (5.8a). 
Its general solution is therefore 


D (t, T) = Ao(r) cost + Bod r) sin t. (5.25) 
The IC’s in (5.24a) yield 
Ao(0) =0, Bo(0) = 1. (5.26) 


Note that we get no information on Ao(T) and Bo (r) from the lowest order 
problem except their initial values. Inserting (5.25) into the PDE in (5.24b), 


we have 
Dit + v= 204 + 40) sin t — (Bh + B) cost). (5.27) 


The Key Argument. The right side of (5.27) gives rise to particular 


solutions for 5 of the form 
Ci(r)t sint and Di (r)t cos t. (5.28) 


To suppress such secular terms, we must set the coefficients of sin t and 
cost to zero in (5.27). Thus Ap and Bo must satisfy the ODE’s 


Ag ＋ A0 = 0, B50 ＋ Bo = 0. (5.29) 
The solution of these equations, satisfying (5.26), are 
Ao(r) O, Bo(r) Se“. (5.30) 


We have now not only guaranteed that b, as a function of t, is bounded, 
but have also determined b explicitly. Thus. 


= +0(€) 
Se sint Ole) 
Se sint + Ode), (5.31) 


for all t >0. 
Exercise 5.4. Compute ò (t, T) explicitly. 


Exercise 5.5. Compare (5.31), with the ù term included, 
against the exact solution given = Exercise 5.1 and explain the 


significance of the term re in v. 
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Exercise 5.6. Compute the energy E(t, €) exactly, using the re- 
sults of Exercise 5.1. Compare your expression with the heuris- 
tically derived approximation (5.18). 


Two Scales and a Strained Variable Combined. The exact solu- 
tion given in Exercise 5.1 to the IVP (5.4) shows that the origin of the slow 
time 7 = et is the factor e. The exact solution also reveals something 
that we failed to anticipate: the fast time is not quite t but rather the 


strained time 
T=tV1-e?. (5.32) 


Damping alters the period of oscillation. Failure to introduce a strained 
time results in the appearance of terms such as er“ in v. (See Exercise 
5.5). These do no real harm. Nevertheless, since (rue) = nen = 
endn -i), the max |v (t,7)| grows with n. Moreover, since the real part of 
rneir is r"cost, e does have the formal appearance of a resonance. 
Thus, for aesthetic reasons if no others, it would be nice to suppress the 
T"e-7’s, You are asked to do this in Exercise 5.7 which represents a good 
test of your mastery of Chapter IV and what we have covered so far in this 
Chapter. 

Exercise 5.7 will illustrate another important principle in perturbation 
theory. 


There is no unique way to distribute the original independent 
variable t between the the two new variables T and 7. The rule 
of thumb is to use any arbitrariness to simplify the expansions 
as much as possible—an admittedly subjective criterion. 


Although the exact solution to the IVP (5.4) suggests (5.32) as the simplest 
fast scale, Exercise 5.7 will not automatically produce this result. Rather, 
if it is assumed that 

TS (IT Miet, (5.33) 


then it will be found that there is no way to determine A1. The explanation 
is that the exact solution can also be expressed in terms of a diff erent fast 


time, say, 
T, = tVI- e - iet tVI - e — Mr, (5.34) 
in which case (5.5) takes the form 


VI - euT., 7) er sin(T, + Air) (5.35) 
= e~*" [cos(A1 7) sin T, + sin(Àı7) cos T,]. 
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Thus an arbitrary value of À; is to be expected in Exercise 5.7, but A; = 0 
will give the simplest expression for vo. 


Exercise 5.7. 


(a) In (5.4) set 
r=, T= Ale (5.36) 
and assume that v = v(T,7,€) to derive a new IVP 
analogous to (5.22). 

(b) Assume regular expansions for v(T,7,¢€) and Me) and 
substitute these into the IVP obtained in (1) to get a 
sequence of IVP’s. Carry the analysis far enough to 
suppress the reg term in ò. 

(c) Check your results against the exact solution for 5 
given in Exercise 5.1. 


Nonlinear Damping. We now examine a problem for which there 
exists no exact solution—a cubically damped oscillator (This example is 
Cole’s): 

b ＋ c +v=0, 0<t, v(0)=0, 500) 1. (5.37) 


Before we attempt to solve this equation, we note that the equation and 
boundary conditions are unchanged if we replace v with —v, t with -t and 
e with —e. Since it is known that the solution (5.37) is unique, then this 
solution must exhibit the symmetry 
v(t, e) = —v(—t, —e) (5.38) 

An attempt to solve (5.37) by assuming a regular perturbation expansion 

v(t, e) = volt) + evi(t) +- (5.39) 
leads to resonant terms in vi (t). 


Exercise 5.8. Show that a regular perturbation expansion of 
the solution of (5.37) leads to secular terms. 


An examination of the (approximate) rate of decay of the energy again 
suggests that we assume 


v=v(t,7,¢<, 1 = d. (5.40) 
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With the aid of (5.20) and (5.21), the IVP (5.37) takes the form 
vrt tu + elu? + Quer) + €2(3u20, + vrr) + 3v? + e =0 (5.41) 
v(0,0,€) = 0, (0,0, e) + €v,(0, 0, €) = 1. 


Exercise 5.9. Verify (5.41). 


As with linear damping assume that 
vlt, 7, e) =v (t, 7) tev (, 7) 4 (5.42) 
This series, substituted into (5.41), implies that 


e: vn + - 0, 9 (0, 0) = o, % (0,0) =1. (5.43a) 


3 
er: be +040; 12 94 0, 0 (0,0) =0, ös (0,0)+ L, (0,0) =0, (5.43b) 
etc. 


Exercise 5.10. Find the next equation in the sequence of 
(5.43). 


The solution of (5.43a) is 
v= Ao(r) cost + Bo(r) sin t, (5.44) 


where 
Ao(0)= 0, Bo(0) = 1. (5.45) 


By the symmetry of (5.38), 
v (t, 7) = v (-t, 7). (5.46) 


This symmetry implies Ap = 0. Substituting (5.44) (with Ao = 0) into the 
PDE in (5.43b), we have 


ber + d= —2B} (7) cost — [Bo(r) cos thb. (5.47) 


Then, with the aid of the trigonometric identity (4.28), we may cast (5.47) 
into the form 


Ün + v= —(2B + 30/4) cost — (B0/ 4) cos 3t. (5.48) 
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Clearly, to avoid a resonant term in b, we must set 
250 + 3B3/4 = 0. (5.49) 


The solution of this nonlinear ODE, satisfying the IC Bo(0) = 1, is found 


easily to be i 


Bo(r) = JOA tE (5.50) 
Hence, from (5.44), 
a sin t sint (5.51) 


Mr MFI 
Note that 
1. Damping is algebraic rather than exponential. 


2. The influence of e has been brought into the first approximation so- 
lution v (t r) through the appearance of the slow time 7 = et. 


3. As with linear damping, it was necessary to consider the IVP for 5 
before d could be pinned down completely. With Bo(r) in hand, we 
can now solve (5.48): 


cos 3t 


32((3/4)7 + 1)" (5.52) 


ù (t, T) = Aı(T) cost + 


Exercise 5.11. Carry the analysis far enough to determine ù 
explicitly. Are there any indications that a fast time T' = A(e)t 
should have been introduced? 


Exercise 5.12. Consider the following IVP for a quadratically 
damped oscillator: 


b+elfb+v=0, O<t, (0) = 1, 500) 0. 


Assuming that e is small and positive, look for a solution of the 
form 1 j 

vlt, 7, €) =v (t, 1) + €v (t, T)+ >, 
where 7, the slow time, is to be suitably related to t. Carry your 
work far enough to determine v explicitly. Hint: See Exercise 
4.4. 
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Exercise 5.13. Consider (1.90) page 24, the dimensionless 
equation for a beam under tension on an elastic foundation. 
Set f(z) = 0 and assume that the beam is semi-infinite and 
subject to the BC’s 


y(0) = 1, y/(0) =0, y(z), y'(z)— 0 as z > ov. 


(a) Solve the BVP exactly for O < e < 1/ M. If O <€ <1, 
the solution exhibits a very short scale and a very long 
one. What are these? 

Introduce a short length X and a long length € to 
account for the two disparate length scales found in 
(a) and assume that y = AX, &, e). Use the two-scale 
method to determine a uniformly valid first approxi- 


0 
mation 9 (X, &) to 9, e). Compare your result with 
the exact solution. 


(b 


— 


Chapter 6 


The WKB 
Approximation 


The WKB?! method is used to approximate the solution of differential equa- 
tions of the form 


y" +u?r(z)y = 0 


where w? >> 1. Before explaining the method, we develop an example of 
the how such an equation may arise. 

A tightly stretched string of length L under a tension T is shown 
in Fig. 6.1. The string has a mass/length (C), where E is distance along 
the string from the left end. If a particle initially at € undergoes a small 
transverse displacement uE, t), where t is time, then its equation of motion 
is 

w Fw 


A derivation of (6.1) may be found in a number of texts. 


Exercise 6.1. Reproduce such a derivation. 


1The initials stand for Wentzel, Kramers, Brillouin. Jefferys also discovered the 
method in the 1920's to deal with quantum mechanical problems. For this reason the 
method is sometimes called the “WKBJ” method; however, according to Steele (“Ap- 
plication of the WKB Method in Solid Mechanics", Mechanics Today, Vol. 3, Ed. S. 
Nemat-Nasser, Pergamon Press, 1976), priority should go to Carlini (1817), Green (1837) 
and Liouville (1837), but the name WKB is now traditional. 
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r(x) 
w(, t) 


ta 


Figure 6.1: A tightly stretched string of variable density. 


Standing waves are, by definition, motions of the string such that 
each particle moves harmonically in time, i.e., 


w(E, t) = H (e) cos Mt, (6.2) 
where Q is the frequency of vibration. Substituting (6.2) into (6.1), we 
obtain the ODE ow 


— 2 = U. b 
Ter Y =0 (6.3) 


We shall require that the ends of the string remained fixed. Thus the 
solutions of (6.3) are subject to the BC’s 


W(0) = W(L) =0. (6.4) 
Nondimensionalization. Let 
€=Lz, W=Ly, p pr, w’ = pL" /T, (6.5) 
where p = max p(é), 0 < Ẹ < L. Then the DE (6.3) and BC’s (6.4) reduce 
to 
Tor) y = O, O< z<1, y(0)= 1) =0. (6.6) 


The Eigenvalue Problem. The BVP (6.6) has the trivial solution 
y(z) = 0. But if r(z) > 0 on (0,1), then there exist an infinite number of 
discrete values of w, 0 < wo < w; < ..., such that for each non-negative 
integer n, (6.6) has a non-trivial solution n (2). wn is called the n“ nat- 
ural frequency of vibration and n (2), the associated mode.” Finding the 


2See, for example, Burkill, The Theory of Ordinary Differential Equations, 2nd Ed., 
Oliver and Boyd, 1962, Chapter IIL 
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Figure 6.2: Qualitative behavior of the eigenfunctions of Eq. (6.6). 


pairs {wn, Yn} is called an eigenvalue problem (EVP). Values of À = w? for 
which (6.6) has non-trivial solutions are called eigenvalues and the asso- 
ciated solutions y,(z), eigenfunctions. (Because EVP’s are homogeneous, 
eigenfunctions are determined only to within an arbitrary multiplicative 
factor.) It may be shown, as indicated in Fig. 6.2, that y,(z) crosses the 
interval (O, I) of the x-axis exactly n times, i.e., the nth mode has n nodes. 


The WKB approximation for High Frequencies. Closed form 
solutions of (6.6) exist only for the simplest functions r(z). However, if 
wn > 1, then between any two successive zeros of y,(z), r(z) is nearly 
constant (see Fig. 6.2). This observation may be exploited as follows. 

If r were a constant, then (6.6) would have solutions of the form 


ys(z,Q) = Aer C7. (6.7) 


This suggests that we assume 


y(z,w) = e). (6.8) 

Then 
y = wg! (6.9) 
y" = (97 +wg")e9, (6.10) 


and the DE in (6.6), upon division by w?, reduces to 


wg" +97 4+r=0. (6.11) 
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As (6.11) is first order in g’, we set 
g= / hdz, (6.12) 
so obtaining 
wh’ +h? +r=0. (6.13) 
Always, we start with the naive assumption that when a problem con- 
tains a small parameter (“ in this case), its solution is regular in this 
parameter. Thus, we set 
h(z,w) = hol) Th (6.14) 


Substituting (6.14) into (6.13) and equating to zero coefficients of successive 
powers of , we obtain the sequence of equations 


he T T= 0 
2hohı + ho =0, (6.15) 
etc. Since r(z) > 0, the first of (6.15) has the solution 
ho(z) = +ir'/?(z), (6.16) 
and the second of (6.15) yields 
hi(z) = —(1/4)[Inr(2)f (6.17) 
So far, 
y(z,w) = exp[w J hodz + $ k dr + O(Ww~')] (6.18) 


=r7/4(2)explu(+i fe ?(t)dt + O(w™?))], 


where the lower limit of integration a is to be chosen for convenience and 
multiplicative constants of integration have been set to 1. In terms of real 
quantities, the general homogeneous solution is therefore of the form 


y(z,w) r (2) {A cosul | / dt + O(w~?)) 
= 
+ B sinu[ ji r!/2(t)dt + 000 (6.19) 
where A and B are arbitrary real constants. It may be proved rigorously 
that (6.19) is an asymptotic expansions. (See Erdélyi, Asymptotic Expan- 


sion, Dover, 1953 p.79.) When the O-terms are omitted, (6.19) is called 
the WKB approzimation. 
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Exercise 6.2. Show that 


babe) 4 7175 THO — . 51 (6.20) 


Exercise 6.3. Since exp[O(w~")] = 1+O(w~"), the expansions 
(6,19) can be written, alternatively, in the form 
y(z,w) = A(z,w) exp[wg(z)] 
= [Ao(2) )- expli | ra 


Derive the DE’s satisfied by Ap and A; and note the relation of 
A, to hz. 


First Order Solution to the Eigenvalue Problem. Setting a = 0 
and substituting (6.19) into the BC’s in (6.6), we find, in the limit as 
w oo, that 


A=0, Bsinſo { ' r!/2(t)dt] = 0. (6.21) 
0 


To avoid the trivial solution, the coefficient of B must vanish. This yields 
the asymptotic formula 


ee C = 
Wn Wn = Or n co. (6.22) 
Hence, 2 
Yn(2,Wn) ~ Bara / () sinkon | r1/2(4)d¢], (6.23) 
0 


where the B,, are arbitrary constants. 
Comparison with some Exact Solutions. If r(z) =z”, m > —2, 
the exact solution of (6.6) may be expressed in terms of Bessel functions: 


vn (æ, Wn) = Bur / (J), (6.24) 
where 1 
P = Lea (6.25) 


ur is the nth zero of Jp(z), and 


wn = p. (6.26) 
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(See, for example, Hildebrand, Advanced Calculus for Applications, 2nd 
Ed., Prentice-Hall, 1976, p. 153.) From (6.22), 


Ôn = (1/2)(m+ 2)(n + )x. (6.27) 


If m = 0, (6.23) and (6.27) are exact, whereas if m = +1, pl? may 
be found in published tables. ( E.g., Watson, Theory of Bessel Functions, 
2nd Ed., Cambridge University Press, 1962, Table VII). Table 6.1 compares 
values of wn from (6.26) with those given by (6.27). 


Table 6.1 Exact and approximate natural frequencies of (6.6) 
with r 241. 


Exercise 6.4. Show, formally, that with the change of depen- 
dent variable w = vy and a proper choice of the function v, the 


DE 
a(z)w” + b(z)w' + [c(z)+w7d(z)]w=0, 0<z<1 (6.28) 
reduces to the so-called normal form 


y” + [9(z) +w?r(z)Jy=0, 0<z<1. (6.29) 


Exercise 6.5. Use the results of Exercise 6.4 to reduce the 
following DE’s to normal form. 

(a) w” + 2c +w?w =0 

(b) (1-7) 22 +w?w =0 

(c) w” + zw +(z?+u?)w=0. 


Exercise 6.6. If r(z) > 0 on (0,1), show, formally, that the 
general solution of (6.29) is of the same form as (6.19). 
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Exercise 6.7. If e is small and positive, construct a uniformly 
valid first approximation to the linear oscillator with damping 
proportional to time: 


t+edv+vu=0, 0<t, v(0)= 0, 0(0)=1. 


To start, put the differential equation in the form of (6.29), then 
try the naive expansion 


v(t, e) = volt) + ev, (t) +- 


Exercise 6.8. Consider the natural frequencies of a tightly 
stretched string composed of two strings of equal length, welded 
end-to-end, where the mass/length of one of the pieces is four 
times that of the other. It follows that r(z) = 1/4 on the interval 
(0, 1/2) and that r(z) = 1 on the interval (1/2, 1). 


(a) Compute the (dimensionless) natural frequencies by 
solving the DE in (6.6) on the two domains (0, 1/2) 
and (1/2, 1), subject to the BC’s y(0) = y(1) = 0 and 
the junction conditions that y and / be continuous at 
z= 1/2. 

(b) Compute ûn from (6.22) and compare your answers 
with those obtained in the part above. To use a phrase 
of our youth, the WKB approrimation is not “phased” 
by discontinuities in the differential equation! (This is 
a drawback as well, for if we study traveling waves, 
the WKB approximation, unless modified, does not 
predict the reflection that must occur at a discontinu- 
ity.) 


Exercise 6.9. Find a change of independent variable € = (z), 
in terms of integrals of a(z)/b(z), such that (6.28) reduces to 
the alternate normal form 

ü + lo) Ye =0, ù = dw/dé. (6.30) 


Exercise 6.10. Reduce the DE’s in Exercise 6.5 to the form 
(6.30). Compare the advantages (if any) of this form to (6.29). 
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Exercise 6.11 
(a) Show that (6.6), (6.19), and (6.22) imply that 


Wn = On + O(n"). (6.31) 


(b) Use the results of Exercise 6.2 to compute an O(n-1) 
correction to ûn. 

(c) Explain what goes wrong when you attempt to use 
the results obtained in (b) to compute an “improved” 
approximation to wn for r(z) = z or N) = sin z. (See 
Exercise 7.5) 

The WKB approximation for Non- Oscillatory Solutions. With 
the change of parameter 


w= e-, (6.32) 
the DE in (6.6), upon multiplication by e2, takes the form 
ey’ —r(z)y=0, 0<z<1, r(z)>0. (6.33) 


Exercise 6.12. Show that any solution of this DE can cross the 
z axis on the interval [O, I] at most once. Hint: Recall Rolle’s 
Theorem (a special case of the mean value theorem) which states 
that if (21) = y(z2),21 < z2, and if y(z) is differentiable on 
(z1, z2), then there exists an 2. E (z1, z2) such that y (z.) = 0. 
Write e?y” = ry and integrate between appropriate limits to get 
a contradiction. 


The WKB approximation to the two solutions of (6.33) follows imme- 
diately from (6.19) and (6.32): 


ys(z, €) r ( Y) exp[te! J ‘ r!/2(t) dt]. (6.34) 
0 


If 0 < e < 1, then e! is large and the exponential factor in (6. 34) 
associated with y_ decays rapidly from 1 towards 0 as z increases from 
0. See Fig. 6.3. In the WKB approximation to y+, it is convenient to 
have an exponential term that displays the same behavior with respect 
to the right end of the domain (0,1). As an arbitrary constant times a 
homogeneous solution is still a homogeneous solution, we multiply the WKB 
approximation to y by exp(—R/e), where 


Hs f * idt. (6.35) 
0 
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n(x) 
1 
x 
0 1 
Figure 6.3: The WKB approximation to the exponential decaying solutions 
of Eq. (6.33). 

Since i = i 

foia- [Cae f at (6.36) 

0 0 z 


it follows that, in place of y4, we may write 


1 
Ü+ bee r1/2(¢) dt). (6.37) 


z 


The exponential factor now decays rapidly from 1 to 0 as x decreases from 
1. See Fig. 6.3. The WKB approximation to the general solution of (6.33) 
is therefore 

y ~ Ay- + Bis, (6.38) 
where A and B are arbitrary constants. In y_ and gy, we have another 
example of boundary layers—functions that decay rapidly to zero as we 
move from the boundary of a domain into its interior. 

A further simplification of the WKB approximation to decay- 
ing solutions is possible because the behavior of exp( {> r!/?(t)dt) and 
exp( J; r!/2(t)dt) is important only near z = 0 and z = 1, respectively. 
This is clear from Fig. 6.3: away from the end-points of the domain (0, 1), 
the exponentials are essentially zero. 

To allow for the possibility that r(x) is singular (but integrable) at the 
end-points, let us assume that 


ao + O(z)], -2, 0<2<1/2 
* { (1 505 2 a i >-2, 1/2 pf <1 (6.39) 


Then > 
| r!/ 2(t)dt = az?(1+ O(z2)] (6.40) 
0 


i r!/?(t)dt = A(1 - z)*[1 + O(1 — z)), (6.41) 
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where 


p= mee >0, a= Vao/p (6.42) 
a= "2? 50, p= Vov/a. (6.43) 


It follows that, asymptotically, 
Ar exp[—az? /e] BCI A Rp EHC 2)5/ e], (6.44) 


where the constants ag 1/4 and 50 * coming from (6.39) have been absorbed 


into the arbitrary constants A and B. 
The DE (6.33) can not arise in an EVP because its solutions do not 


oscillate. Typically, the auxiliary conditions associated with (6.33) are non- 
homogeneous, for example BC’s such as 


„lim. 204% = a, lim (1 ) =b. (6.45) 


Substituting (6.44) into (6.45), we obtain 
A~a, Bb, (6.46) 
since exp(—a/e) ~ 0 and exp(—ĝ/€) ~ O as € + 0. 


Exercise 6.13. Prove this last statement. Hint: Use l’Hospital 
rule to show that lim. e — 0 for every p. 


9 
8 


Exercise 6.14. Explain why two functions f (e) and g(e) may 
have identical regular expansions but still differ by terms such as 
exp(—a/e). (Thus there is a limit to the amount of information 
that can be extracted from a regular expansion as we saw in the 
Stieljes integral (1.69 ).) 


Exercise 6.15. Use (6.44) to obtain an asymptotic solution to 
the BVP 
ey’ —2(1—z)y=0, 0<2<1, 
; (a ; I 

jim, 2 y= 2, „im (1 19774 
Make a careful graph of the solutions for e = 0.1 and e = 0.2. 
Note that the approximate solutions are singular at z = 0 and 
z = 1. In contrast, the exact solutions must be regular at the 
ends of the interval. 


Chapter 7 


Transition Point 
Problems and Langer’s 
Method of Uniform 
Approximation 


A buckled drill string ! of length L is sketched in Fig. 7.la. Fig. 7.1b 
indicates the forces and moment on a typical cross section lying a distance 
s from the lower end of the string. It is assumed that the loads have rotated 
the section through an angle ¢(s). Let w(s) denote the weight / length of the 
string and V,, the force at its upper end. The vertical force at any section 
is then 


L 
7% = J w(t)dt. (7.1) 
s 
From force and moment equilibrium and strength of materials, we have 
V sin G = dM/ds = Eld?¢/ds?, (7.2) 


where E is Young’s modulus (an elastic parameter) and J is the moment of 
inertia of the cross section about a line through the centroid, perpendicular 
to the plane in which bending takes place. The ends of the string are 
assumed to be free to rotate; hence 


M = Eld¢/ds = O at s =0. (7.3) 
1 Drill string” is the name of pipe used to drill oil wells. These pipe are several inches 


in diameter and thousands of feet long (screwed together in 20 foot sections). As the 
drill string can support a transverse force as well as a tension, it is actually a beam. 
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V, 
4 Vis) 
jj os) 
[ | 7 
(a) @ 


Figure 7.1: (a) A drill string as it appears in the well casing (b) Loads on 
a section of the string. 


Linearization and Nondimensionalization. To simplify the math- 
ematics we assume that 


sing xó =y (7.4) 
and nondimensionalize as follows: 
s= Lz, f =T/W, ê = EI/[L’W, (7.5) 
where Š 
1 w(t)dt (7.6) 
0 


is the total weight of the string. The DE (7.2) and BC’s (7.3) now take the 
form 
ey" — f(z)y=0, O C41, y¥(0)=y(1)=0. (7.7) 
The drill string is under tension as it is lowered. However, after its end 
has hit bottom, the lower part of the string will be in compression while its 
upper part will remain in tension so long as % > 0. The dimensionless ten- 
sion f(x) will have the general shape sketched in Fig. 7.2. It is convenient 
to shift the independent variable so that f(0) = 0. We shall therefore study 
the equation 
627) — f(z)y=0, -a<z <b, (7.8) 
subject to certain BC’s at z = —a and z = b which are of no immediate 
concern. 
If z < —ô < 0, f(z) < 0, so, from (6.19), the solution of (7.8) is of the 


form 
You ~ ()J i cos[e~* F(z)] + C2sinfe~'F(z)]}, (7.9) 


A First Look At Perturbation Theory 83 


Nx) 


Figure 7.2: Dimensionless tension in a drill string. 


Figure 7.3: Sketch of f(z) and the solution of (7.8). 


where = 
F(z) = J |F(&)|1 dt. (7.10) 


If O < ô < z, f(z) > 0, so, by (6.34), the solution of (7.8) is of the form 
Ymon ~ |f(2)| 71 [Cs exp[—e7 1 F(z)] CA exp[e1 F(z)]}. (7.11) 


However, in a neighborhood of z = 0 neither solution can be valid because 
limo (z) / — oo [even though z = 0 is not a singular point of 
(7.8).] Nonetheless, as the solution of (7.8) must be continuous on (—a, b), 
we can infer Fig. 7.3. 

The connection problem is to express Ca and Ca in (7.11) as a linear 
combination of Ci and C} in (7.9) (or vice-versa). This problem must be 
solved before we can solve a BVP. There are two major ways to approach 
the connection problem, typified by methods proposed by Gans and Jefferys 
and by Langer. Gans and Jefferys find an approximate equation which they 
solve exactly, whereas Langer introduces a change of variable and solves the 
resulting exact equation approximately. 


84 Transition Point Problems, Langer Method 


Figure 7.4: Graphs of Ai and Bi. 


The Gans-Jefferys’ Method (1923). Near z = 0, (7.8) looks like 
eY" — 2 (OY = o, (7.12) 


provided that f’(0) # 0, as we shall assume henceforth. This is a Bessel 
equation. We expect (hope!) that in some way, as € — 0, 


Y (2) ~ Ymon, 2 >0 (7.13) 


Y(z)~ Yor, 2 0 (7.14) 


The parameter € as well as the constant f (O) may be scaled out of 
(7.12) by setting 


z = [f (0] 252. (7.15) 
This reduces (7.12) to Airy’s equation, 
d*y 
eg -zY =0. (7.16) 


Facts About the General Solution of Airy’s Equation. As z = oo 
is the only singular point of (7.16), this DE has two linearly independent 
solutions analytic everywhere. These are usually denoted by Ai(z) and 
Bi(z), so that the general solution of (7.16) takes the form 


Y(z) = CsAi(z) + CeBi(z) . (7.17) 


Properties and values of Ai(z) and Bi(z) are given on pages 446-452 and 475- 
478 of the Handbook of Mathematical Functions, Abramowitz and Stegun, 
editors, U.S. Government Printing Office, 1964. 

Graphs of these functions are given in Fig. 7.4. 
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As z — o0, 
4 ee (7.18) 

Bi(z) ~ ob (7.19) 

Ai(—z) ~ -yray in (G+ z) (7.20) 

Bi(-z) ~ = 25 74008 Ge" + z) i (7.21) 


Matching via an Intermediate Variable. To apply Jefferys’ method, 
we express both z and z in terms of an intermediate variable ij such that if 
we hold n fixed and let € — 0, then z — œo while z — 0. To achieve this let 


z=Ke%n, z= Ên, K, a, f > O, (7.22) 


where K a, and 2 are constants to be determined. Substituting (7.22) into 
(7.15), we get 


Pn = KITO) l etn. (7.23) 
If (7.23) is to hold for ņ fixed and all €, we must set 
-a+f=2/3, K LH. (7.24) 


Any choice of q and £ satisfying (7.24) will do; for symmetry, we set - = 
8 = 1/3. Then from (7.22), 


se [POPC ey, aay. (7.25) 


We first match (7.11) and (7.17). In this case ņ > 0, so that from (7.18), 
(7.19), and (7.25), 


Y (2) I/ , 
Cs Ail f()) 1/5 + CsBil. . |] (7.26) 
(02/3) 7760072 ˙¹—7(e 55 / 
EE O A U A 5 c. 
while from (7.11) and (7.25), 

Ymon(z) = Imon(e / 

= Cs|f (èn) s F (e"/n)] (7.27) 
+C4| f(n) expl F (n). 
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Now by Taylor’s theorem, 
f(x) Fo) r + Ol). 


Hence 
Fr) = F + Olè?) , 


so that from (7. 10), 
F(z) = 37e $ O(z5/2), z>0 


2 
a qf Ons 


(7.28) 


(7.29) 


(7.30) 


Thus 
—(2/3) f!(0)1/2¢-1/27 3/2 +... 


If (7.26) and (7.31) are to agree, then 


C3 = Cs 
f'(0)1/4¢1/12 aa 271/2 f'(0)!/12e71/12 


or 


10 1/6 
= Eo. 


C3 


Exercise 7.1. Find the relation between Ca and Cs. 


Exercise 7.2. Show that (7.14) implies 


C2 z (27) 72 Cs — Ce 


Exercise 7.3. Obtain asymptotic solutions to the EVP’s 
(a) ey” —[sin(wz/2)]y=0, y(—1) = y(1) = 0. 


(b) ey" —(tanhz)y=0, y'(—1)= y (1)=0. 


(a) -H ( ate} (7.34) 


(7.32) 


(7.33) 
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Langer’s Method. Our treatment of (7.8) for €? < 1 required us 
to construct approximate solutions on the three domains —a < z < —6, 
— < z < 6, and 6 < z < b and then to relate the various constants of 
integration by a matching procedure. The approximate solution on (—6, ô), 
(7.17), was a linear combination of the well-studied Airy functions, Ai(z) 
and Bi(z). 

It was Langer’s idea to attempt to describe the behavior of the solution 
of (7.7) over the entire domain —a < z < 6 in terms of Ai and Bi. His 
method consists of introducing new dependent and independent variables 
such that, except for small terms, (7.8) takes the form of (7.16). Langer’s 
work was initiated in the early 1930’s and completed in the late 1940’s. The 
two-scale method, which first appeared about 10 years later, turns out to 
be an ideal way of arriving at Langer’s results. 

Recall that to study the behavior of (7.8) near the transition point z = 0, 
we approximated f(z) by zf’(0) and then introduced the new independent 


variable 
z = 7/3, f'(0) 32, (7.35) 


Compared to z, 2 is a “fast variable” (in Cole and Kevorkian’s terminology), 
since a change in z of O(e?/%) produces a change in z of O(1). The results 
of the preceding section and the form of (7.35) suggests that we attempt to 
express the exact solution of (7.8) as a function of 


= 77/3 9(2), (7.36) 


and z. To fix g(z), we shall require that, to a first approximation, we obtain 

a DE in terms of the fast variable € alone, identical to the Airy equation. 

Finally, we shall make g(0) = 0 so that the transition point occurs at € = 0. 
Thus with 2 = /s we assume that 


y= M, E, G) (7.37) 
Then 
d 0 
4 g gd . tO y (7.38) 
Y" = Yeo + 287g! yee + B19" ve + H- ye. (7.39) 
Inserting (7.39) into (7.8), and multiplying through by € = gg, we obtain 
90“ yee + B9(2g'Yog L 9“ 16) + Bayes — Ef(z)y = 0. (7.40) 


To make (7.40) take the form of (7.16) as 5 — 0, we set 
99% =f. (7.41) 
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Upon division by (7.41), (7.40) reduces to 


vee + A[(2/9')uee + (9% “) vel + 87 (1/9? vee EY = 0. (7.42) 
Since g has the same sign as f and g(0) = 0, the solution of (7.41) is 
a(2) = sgn(z)((3/2) | t lu (7.43) 


Now assume an expansion for y of the form 


YE, z, 8) =y (E) T y (€,2) +... (7.44) 


Substituting (7.44) into (7.42) and equating to zero the coefficients of suc- 
cessive powers of H, we obtain the sequence of DE’s: 


Yee e Y= 0 (7.45) 
Vee —€ Y= (/ Yer -N) Ye, (7.45b) 


etc., the general solution of (7.45a) is 


Y (€, 2) = Co(z)Ai(é) + Do(z)Bi(é). (7.46) 


To determine the functions Co(z) and Do(xz) we must consider (7.45b). 

For simplicity we determine Co(z) only; Do(z) may be found in a strictly 
analogous fashion. Inserting (7.46) with Do = 0 into the right side of 
(7.45b), we have 


Wee e Y= —[(2/ 9") Ch + (0""/01?)CoIA‘(€). (7.47) 
Observe that if A(€) is a solution of Airy’s DE, then 
(€A)" — €(€A) = 24. (7.48) 


Thus the right side of (7.47) will produce a particular solution of the form 
F(z)€Ai(é) unless the coefficient of Ai’ is zero. To avoid such ”secular” 


terms in 1 (E. 2). we choose Co(z) so that 
(2/91)Cy + G D = 0. (7.49) 
This equation has the solution 


Cole) = aggre /, (7.50) 
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where yo is an arbitrary constant. Hence, from (7.36), (7.43), and (7.50), 
y = olg) / Ai Vg O), (7.51) 


where g(z) is given by (7.43). It may be shown that the error estimate 

is rigorous. (See Section 12.8 of “Asymptotic Methods” by F. Olver in 

Handbook of Applied Mathematics, Ed. C. Pearson, Van Nostrand, 1974). 
As z — +00, we have, using (7.18), 


v~ sere A | oa, (7.52) 


which agrees with (7.10), while as z — co, we have, using (7.20), 


1/12 * 
i e iae | ON d+ aya, (% 
which agrees with (7.9). 
Exercise 7.4. Compute Do(z) in (7.46). 


Exercise 7.5. Carry the computations far enough to determine 


y (2, C). Is there any evidence that a more general fast variable 
of the form € = 6~1g(z, p) should have been introduced? 


Exercise 7.6. Consider the EVP 
y" +w(sinzx/2)y=0, y(0) = y(1) =0. (7.54) 


The WKB method yields (6.22) as a first approximation to 
the n* natural frequency wn, but fails, as you saw in Exercise 
6.10(c), to give a correction. Use Langer’s method to determine 
this improved approximation. 


Chapter 8 


Introduction to Boundary 
Layer Theory 


Boundary layer methods in nonlinear problems were initiated by Prandtl 
near the turn of the century. Prandtl noted that when a fluid of low viscosity 
such as air or water flows about an obstacle, the ratio of viscous to inertial 
forces is small everywhere except in a narrow layer near the boundary of the 
obstacle. Using this observation, he was able to simplify considerably the 
analysis of the governing Navier-Stokes equations. His idea was that there 
is a region far from the obstacle where the flow is essentially the like the flow 
of an inviscid fluid. On the other hand, near the obstacle, where viscosity 
plays an important role in making the velocity equal to zero on the surface, 
the velocity changes much more rapidly along a perpendicular to the surface 
than along the surface itself. This suggests that we introduce the boundary 
layer as ashort interval within which the solution of the differential equation 
changes very rapidly. As we shall see, these intervals of rapid transition may 
occur at the ends of the domain (layer near the boundary), but they may 
also occur in the interior of the domain (a transition point or shock wave). 
Many of the results obtained from using boundary layer ideas can also be 
obtained by using the two-scale method. However, the methods themselves 
are quite different, and you will find that sometimes one method is easier 
to use than the other. 

In Chapter 3 we considered a typical though simple boundary value 
problem, namely, 


ey TY +y=0, 0<2<1, y(0)=0, y(1)=1. (8.1) 
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Figure 8.1: The functions ee~* and -e / and their sum 5̃. 


The approximate solution given there was 
J= efe? e ), (8.2) 


a sum of two functions whose graphs are indicated by the dashed lines in 
Fig. 8.1. The solution of (8.1) exhibits a boundary layer of width O(e) at 
the left end of the interval [O, II. 

The aim of this chapter is to develop ways to obtain, systematically, 
approximate solutions to BVP’s that exhibit boundary layers. To lay the 
footings for later perturbation expansions, we now introduce, via a few 
simple examples, a method of obtaining first-approximation solutions. 

The Boundary Layer Method. The first step is to see what happens 
if we set the small parameter in our problem to zero. Doing this for (8.1), 
we obtain the reduced DE 

Y'+Y =0, (8.3) 


whose general solution is 


Ce-. (8.4) 


As (8.4) contains only one arbitrary constant, C, we cannot satisfy both 
boundary conditions in (8.1). Looking at the graph of 5; in Fig. 8.1, which 
we expect to be a reasonable picture of the solution of (8.1) if e is small, 
we see that the solution changes very rapidly near the left boundary. To 
exploit this observation mathematically, we introduce a new boundary layer 
variable, 


8 /e, (8.5) 
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whereupon the DE in (8.1), when multiplied by €, becomes 


d*y dy 
= +s =0. 8.6 
de d (8.6) 
Setting e = 0, we obtain 
dyst dre 
der de = (8.7) 
whose solution is 
yBL = A+ Be-. (8.8) 


From (8.8), we see that the boundary layer will be on the left because the 
term e goes to zero as we move away from the left boundary. [If there 
were a negative sign between the first two terms in (8.6), the boundary 
layer would have been on the right.] 

We must now try to combine (8.4) and (8.8) to make an approximate 
solution of (8.1). As the boundary layer is on the left, we may find C by 


forcing Y (1) = 1. Thus 
Ne (8.9) 


This is called an outer or interior solution (meaning that it is valid in the 
interior of the region, away from the boundary layer). The solution (8.8) 
must therefore meet the boundary condition on the left, that is, 


yet = A(1 - e$). (8.10) 


The constant A is determined by making (8.9) and (8.10) match at the edge 
of the boundary layer. Heuristically, we want the limit of the outer solution 
as z approaches the boundary layer to equal the limit of the boundary layer 
solution as € approaches the outer region. In symbols, we want 


lim Y = im YBL, (8.11) 
which, by (8.9) and (8.10) yields 


Aze. (8.12) 


In certain problems, matching is more complicated and must be carried 
out in a more formal way by introducing an intermediate variable—call it 
z—similar to that used in Jefferys’ method for the transition point problem. 
To demonstrate the method with (8.9) and (8.10) we set 


r=6°z, =z. (8.13) 
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We want a to be positive and J to be negative so that, for fixed z, z will 
go to zero and & will go to infinity as € goes to zero. Because C = /e, we 
have 

B=a-1. (8.14) 


For symmetry, we take a = 1/2, 6 = 1/2, substitute (8.9) and (8.10) into 
(8.11) and get! 


limee~V = lim A(1 — e~*/V*), (8.15) 
which implies that 
Ae. (8. 16) 
Our approximate solution is now given by two pieces: 
_ f el—e-*/*), z< Ole) 
* { ee z Ole) (8.17) 


Equation (8.17) is awkward because we must decide when to stop using 
one expression and to start using the other. A uniformly valid solution 
may be found by adding the two solutions and subtracting the part which 
is common to both. This common part is simply the common limit given 
by (8.15). Thus 


Yuniform = YBL + Y — common part 
= e(1—e~*/*) + ee™® —e (8.18) 
= ele — . ae 


As a slight modification of our original problem, (8.1), let us consider 
e'-y +y=0, 0<z<1, y(0)=0, y(1)=1. (68.19) 


The first thing we notice is that a boundary layer will form on the right 
because the non-constant solution of cy” — y' = 0 decreases as z decreases. 
The reduced equation, 


-Y'+Y =0, Y(0)=0, (8.20) 
has the solution 
Fan (8.21) 
Introducing the boundary layer coordinate 
=- /, (8.22) 


1Recall that lim is short for “the limit as € + 0”. 
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which is chosen to vanish at the right boundary (z = 1), we obtain the 
boundary layer equation 


dyn dyBL 
— SUBL ii), ii 8.23 
dé? dé yBL( ) ( ) 
with solution 
YBL = 1- B(1 - e£). (8.24) 


If we match (8.24), as we move away from the boundary layer, to (8.21) as 
we move toward the boundary layer, we find that B = 1. Therefore, 


paasa me UE, (8.25) 


Here, the common part, as well as the outer solution, is zero. 
An interesting variation on this type of problem is one which does not 
have a boundary layer at either end. Consider 


ey! +zy' + 229 2 0, 14 241 (8.26) 


„( 1) =a, 1) = B. (8.27) 

There is no boundary layer on the left (x = —1) because the sign between 

the first and second term of the DE is negative, but neither is there a 

boundary layer on the right (z = 1) because the sign there is positive. The 
reduced equation 

zY’'+27Y =0 (8.28) 


has the solution ` 

Y = Ce (8.29) 
We can meet both boundary conditions in (8.27) by choosing C differently 
for z < 0 and z > 0. Thus 


= de le- 1, 1 4 0 
Y= l ge be- l= — (8.30) 


A sketch of (8.30) is shown in Fig. 8. 2 for a and g positive. 

It is clear from Fig. 8.2 that the first and second derivatives of the exact 
solution of (8.26) and (8.27) must be large near z = 0. To study this 
behavior we introduce the new variable 


e = M (8.31) 
in terms of which (8.26) reads 
@? 
Te +e 4 27% 0. (8.32) 
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P x 
-1 0 1 
Figure 8.2: Solution of the reduced DE (8.28) meeting the BC’s (8.27). 


Setting € = 0 produces a boundary layer equation with solution 


00 
yo. 475 I e / d, (8.33) 


[For simplicity let J (C) = fe° e~/?"dt. Then 2J(0) = J(—co) = V2z]. 
Using the quick matching procedure 


— YBL = E. Y 
lim yar = lim Y, (8.34) 
=- 2—-0- 


we obtain equations for A and B: 


A pe!l? 
P = (8.35) 


Solving (8.35) for A and B, we produce the following uniform solution: 


Se] e — B)I(z/Ve)/V2x, z> 0 
Panitorm = | der/ ei- HIV Vr 1], 2 <0 
(8.36) 
At 1 = 0 each expression on the right of (8.36) has the value (1/2)(a + 
be. 

The reader should notice that the first two examples in this chapter were 
ones that we could have solved exactly whereas the last example would have 
been quite difficult to do exactly. A more complicated example, illustrating 
matching in a problem of physical origin, is given in Chapter 9. 

The Two Scale Method. If we examine (8.2), we see two scales; a 
slow one z and a fast one = 2/e. If these are introduced formally into 
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(8.1), as was done in chapter 5, we obtain 


Pp 2 Hy , 1er) , by ley, 
( okt age) a ih en 
y(0,0,e)=0, y(l,e7*, ) = 1. (8.38) 


Multiplying (8.37) by € and setting 


56. . t) =Y (2, Le, +=, (8.39) 
we have for the coefficients of e° and e}, 
0 0 
Vee + Ue 2 Yee . — Y . (8.40b) 
The solution of (8.40a) is 
y= A(z) + B(z)e~®, (8.41) 


which reduces (8.40b) to 
Yee + Y= (B! — Bee —(A' + A). (8.42) 


To suppress terms like ge-“ and E in the particular solution of (8.42), we 
set the coefficient of e-, and €° to zero. This yields 


A'+A=0, B/-B=0. (8.43) 
The boundary conditions (8.38) imply that 
A(0) + B(0)=0, A(1)=1. (8.44) 


Solving (8.43) and (8.44), we obtain from (8.41) 
y= ele eve). (8.45) 


Note that we assumed the boundary layer was on the left when we set 
€ = z/c because € = O(1) in a neighborhood of z = 0. If the boundary 
layer had been on the right, we would have failed. 

The Two-Scale Method for Variable Coefficient DE’s. Consider 
the variable coefficient DE and BC’s 


ey + a(z)y +b(z)y=0, O< e<1, y(0)=0, y(1) =1, (8.46) 
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where a(z) is continuously differentiable and positive on (0,1). If e is 
small, a(z) changes little over an interval of length O(e). Applying the 
same heuristic reasoning as used in Chapter 3, we might expect a uniformly 
valid, first approximation solution of the form 


g ele (. _ ¥(z)), (8.47) 


where ā is some average value of a taken over a small interval near z = 0. 
To obtain more quantitative information, we assume that the solution 
of (8.46) depends on a slow variable z and a fast variable 


€=c7'g(z), where g(0)=0, (8.48) 
and proceed formally: 


y= , &, e) 
y = ye TE (8.49) 
y" = tee + Qyege 1g Teig“ 4 Nec 9% 


Substituting (8.48) and (8.49) into the BVP (8.46), we get, after multiplying 
by e all equations? 


9 ee +ag'ye + e(2yseg' + yeg” + ays + by) + 6 = 0 (8.50) 
y(0,0,e)=0, y(1,e719(1), e) = 1. 
We now assume an expansion for y of the form 
0 1 

Mr, &, e) =y (T, ) Te (r, c) , (8.51) 
substitute this into (8.50), and equate to zero coefficients of successive pow- 

ers of e. This yields the sequence of BVP’s 

0 0 0 0 
9“ Veg Lag = 0 1 (0,00 0 (1,99 =1 (8.524) 

22 oi 0 ia BE mie g 0 
9 Yee Tag“ Ye +2 Vee g + % 9“ +aY: +bY=0 
1 

Y (0,0) 0, y(1,e749(1)) =0, ete. (8.52b) 
2 We could have, by a change of the dependent and/or independent variable, reduced 
the DE in (8.46) to one of the normal forms discussed in Chapter 6. We do not because 


the success of boundary layer methods does not depend on this reduction which may not 
exist for nonlinear or higher order DE's. 
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To make (8.52a) as simple as possible and to meet the initial condition 
in (8.48), we set 


g= J "unit ta Fai (8.53) 
The partial differential equation in (8.52a) then reduces to 
C o 
Yee + Ye=0, (8.54) 
which has the solution 
Y (z,€) = Ao(z) + Bo(z)e~€. (8.55) 
The first BC in (8.52a) yields 
Ao(0) + Bo(0) = 0. (8.56) 


At z = 1, e~€ = exp[—e7? 15 a(t)dt] S Ce- Ke, C, K constant. Thus e“ 
is a transcendentally small term and the second BC in (8.5 2a) reduces to 
40 (1) = 1. (8.57) 

To determine Ao(z) and Bo(z), we must . the form of the solu- 


tion of (8.52b). With 5 given by (8.53) and y by (8.55), there follows, 
after division by 9“ 72 d 


-ę [Bo (4 — 2B) Aj a) 
sT . eee 
Yee + Ue=e l zi ( . (8.58) 
To avoid “resonant” terms of the form F(z)€e~£ + G) in the solution 


for Y, we must set to zero the coefficients of e~€ and &“ on the the right of 
(8.58). Thus, 


44 + bAo = 0, (8.59) 
which in view of (8.57) has the solution 
1 
A(z) = exp{ T lot) /at lr, (8.60) 
while 
aB + (a! — b)Bo = 0 (8.61) 


has the solution 
Bo(z) = [Co/a(z)] exp d [b(t)/a(t)]dt }. (8.62) 
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To determine the constant Co in (8.62), we apply the BC (8.56). With 
Ao(z) given by (8.60) and 


s(2) = [ Oaa, (8.63) 
we find that 
Co = —a(0) exp[f(1)]. (8.64) 
So far, then, 
y = 4 (€,2)+O(6) 
= explf(1) - T f (8.65) 
-{a(0)/a(z)]expl-e! [ a(t)dt + f(z) + f(1)] + Ole). 


Exercise 8.1. Determine Y (£, z) for the BVP 
ey" T (ITT T2 = O, y(0)=0, 1) =1. 


Exercise 8.2. Determine Y (C, z) for the BVP 
ey +2y'+27y=0, y(0)=0, y(1)=1. 


(a) Try the problem using both the two-scale method and 
the boundary layer method. Which method do you 
prefer? 

(b) Show that the differential equation can be reduces to 
the normal form 


Pw as 
dz? —(a+ 72 jw =0 


where a = y1 — 4% (parabolic cylinder function). 


Exercise 8.3. The BVP 
e +a(z)y =0, y(0)=0, y(1) =1, 


can be solved exactly, since the DE is first order in 7. Compare 
the solution with (8.65). 
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Exercise 8.4. Find the expression analogous to (8.65) for the 
BVP 


ey” Ta) / T b) Y =, ¥(0)=1, (1) 0, 
where a(z) is differentiable and positive on (0, 1). 


Exercise 8.5. Consider the BVP 


Ley = e dr) T= O, 0<2<1, 
500) 1, ¥(0)=0, y(1)=y(1)=0. 


Suppose that a(z) > 0 on (0, I) and that two independent sol u- 
tions of Loy = 0 are known, call them U(z) and V(z). Find a 
uniformly valid approximation to the solution using the bound- 


ary layer method. 


Exercise 8.6. Suppose that a beam under tension rests on a 
foundation with a modulus that varies linearly along the length 
of the beam. Nondimensionalizing as in Exercise 1.14, assum- 
ing no external distributed load, and considering a semi-infinite 
beam under the same BC's as in Exercise 5.14, we arrive at the 
BVP 


[y+ (1+ mr) / = O, O K, Om 
500) = 1, 0) o, Nr), 7) 0 as z O. 


If O < e <1, there are, as in Exercise 5. 14, two disparate scales. 
Introduce appropriate short and long lengths, X and &, and de- 


ter mine a uniformly valid frst- approximation ; (X, O to y(z,€). 


Matched Asymptotic Expansions. The two-scale method can fail 
if closed form solutions cannot be obtained for the associated PDE’s. In 
this case, a solution by a so-called matched asymptotic expansion proce- 
dure may sometimes work. Moreover, it may lead to a uniformly valid 
first approximation to a solution more easily than the two-scale method. 
As applied to (8.46), the method of matched asymptotic rests upon the 
recognition that the solution has distinctly different behavior in different 
subintervals of [O, IJ. In the “inner” or boundary layer region, [O, O(¢)], the 
solution is such that 


ey" + c) + O(c) =0, (8.66) 
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whereas in the “outer” region, [O(e), 1], the solution of (8.46) is such that 

a(z)y + b(z)y + Ole) = 0. (8.67) 

The form of (8.67) suggests that in the outer region, we assume an 
expansion of the form 

y = (, e) = Yo(z) Te) T (8.68) 


Substituting (8.68) into (8.46) and equating to zero the coefficients of like 
powers of e, we obtain the sequence of DE’s, 


a(z)Yo + MY) Vo = 0 (8.69a) 
(r) + GO) U = —-Yel1, k = 1, 2, , (8.69 b) 
which has the sequence of general solutions 
Yo(2) = Coexpl-f(2)}, fe) | G (8.70) 
Vi(z) = {a +Co [3 (2) - a at) ewir (8.70b) 
etc., 
Application of the BC y(1) = 1 yields 
Co = exp(f(1)] (8.71a) 
1 2 ! 
c= | al (0 = (8.71b) 
etc. 


In the inner region, we set 
z=, (8.72) 
so that € = O(1), and replace the variable coefficients in (8.46) by their 
Taylor series Pia This puts the DE into the form 


Ey (ao + areg + aeg? + 0K 
+ €(bo + biet + bzeze: T. ) =0. (8.73) 
Assuming an expansion of the form 


y=Y(E,¢) = YE) + F(E) + (8.74) 
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in the inner region, substituting this into (8.73) and equating to zero coef- 
ficients of like powers of e, we obtain the sequence of DE’s 


Vp" + aÑ% =0 (8.75a) 
Yj’ + aofi = -aÝ — boo, (8.75b) 


etc., which has the sequence of solutions 
Ñ (€) = Bo + Coe% (8.76a) 
Ñi (€) = Bı + Cie — (Bobo/ao)€ 
+ Co(bo/ ao — a1 /ao — 1/2a£)€e" e, (8.76b) 
etc. Applying the BC y(0) = 0, we obtain 
Bo+tCo=0 (8.77a) 
Bi ＋ Ci = O, (8.77b) 


etc. 

The Matching Procedure. To obtain additional conditions for the 
three sets of constants A+, Bi and Cx, we must match the inner and outer 
expansions. To this end, we choose an intermediate variable 7 such that 
if 7 is fixed and e — 0, then E — oo and z — 0. A convenient choice, 
consistent with (8.73), is 


e Sin, z= Bn, where g 6. (8.78) 
We now express both the outer and inner solutions in terms of 7: 
Y = Yo(6n) + 67Yi(8n) + --- 


‘= ef) — (bo /a0) Bnet) + O(8?). (8.79) 
Ý = Ý( ½% NC K 
= Bo — Bo(bo/a9)An + O55). (8.80) 


Comparing (8.80) with (8.80), we see that Bo = ed). To determine Bi, it 
is necessary to consider the O(4?)-terms. 
The composite expansion is defined to be 


v. = Y +Ý common part. (8.81) 
In the present case, the common part, to lowest order, is e/(), Hence 


Y, ef- _ ¢~902/¢+4(1) 4 O(e), (8.82) 
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This is essentially the same result we obtained by the two-scale method. 
See (8.65). 

Corner or Interior boundary Layers. In the preceding analysis, we 
assumed that a(z) > 0 on (0,1). If this is not the case, e.g., if a(z) has 
a zero on (0,1), interior boundary layers may exist. The following simple 
problem illustrates such a phenomenon. 


ey! + 227 —2y=0, -1<2<1, y(-1)=A, y(1)=B. (8.83) 


[The interval (-1,1) has been chosen rather than (0,1) to simplify the sub- 
sequent algebra.] A study of the solution of (8.83) is useful in more general 
problems, just as the study of Airy’s DE proved to be useful in the study 
of general transition point problems. See Exercise 8.8. 

As usual, we first consider the reduced DE 


-Y=0, Ic I, (8.84) 


which has the solution 
Y = Ce: (8.85) 


Examining (8.83) for boundary layers, we note that there can be no 
rapidly decaying solution near the left- or right-hand boundaries, since 2r— 
the coefficient of / — is, respectively, negative or positive in these regions. 
Yet if ey” were small everywhere, we would have y ~ Y, which does not 
contain enough arbitrary constants to satisfy the BC’s. As (8.84) has a 
singularity at z = 0, we would be justified in taking a different value for C 
in (8.85) for z < 0 than for z > 0. These values of C can be chosen so that 
both boundary conditions in (8.83) are met. Nevertheless, the resulting 
function is not a solution of (8.83) because its first derivative is different on 
either side of z = 0. 

In the present example, the situation can be analyzed exactly by first 
noting that (8.85) is, in fact, an exact solution of (8.83). This reduces 
the determination of the complete solution of (8.83) to quadratures, via a 
reduction of order. 

To simplify algebra, let 


4 K. (8.86) 
Then the DE in (8.83) reads 
d) „pdy — “il 
de: N — 2 0, O HCK (8.87) 


Note that the scaling (8.86), while removing € completely from the DE 
(8.87), has introduced ¢ into the solution domain. This is an illustration of 
the third point in the introduction to Chapter 3. 
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Since y = € is a solution of (8.87), we set 


= Eule). (8.88) 
In terms of u, (8.87) reduces to 
eu“ + 2(1+ E) u“ = 0. (8.89) 
This is a first order linear DE in u. Its solution is found easily to be 
ul = cget, (8.90) 
and so 
imam [ at 
Sei- ei Ele- — 2 fi r e-t dt) (8.91) 


= ch Teile le- Vr + erf(€)], 


where erf(£) is the well-studied and tabulated error function. Therefore, 
by (8.88), 


y = czẽ Heile Vr + £ ert 
= Agr + Aille Vr e + zerf(z/e)}. (8.92) 


Imposing the BC’s in (8.83), we have 
A= Az — A; + transcendentally small terms (TST) in e 


B = Á A1 + TST ine. (8.93) 
Thus, 
— [5 — Bee tials o|- 34 ce J + TST. (8.94) 


The sketch of (8.94) in Fig. 8.3 shows the existence of an interior or “corner” 
layer of width O(e) centered at z = 0. 
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Figure 8.3: The solution of (8.83) which exhibits a corner layer. 


Exercise 8.7. Discuss the solution of 


ey" + ary —y=0, 12 4,1, y(-1)=A, 1) = B. 


Exercise 8.8. Discuss the solution of 
ey" + f(z)! —y=0, -l<2<1, x-1)=A, (1) B. 
where f(0) = 0 and zf(z) >0, 2 K 0. 


Chapter 9 


Cables and Cells: Ancient 
and Modern Problems 


In this chapter, we demonstrate the solution of two problems which occur 
in nature. Though Prandtl developed boundary layer methods for solving 
fluid mechanics problems, these problems are a little too elaborate for this 
book,! but Prandtl’s ideas may be illustrated in other ways. To hint at the 
efficacy of boundary layer techniques for solving nonlinear problems, we 
have chosen two problems, one old, one new, both of practical importance, 
we think, and both without exact, closed-form solutions. 

The Shape of a Hanging Cable with Small Bending Stiffness. 
This problem springs from the much older Catenary problem: to find the 
shape of a hanging chain, i.e., a cable with no bending stiffness. This might 
well be called THE PROBLEM of late 17th Century physics. Though now 
stated and solved in terms of the hyperbolic cosine in nearly every calculus 
book, the problem challenged the greatest mathematicians of the time. 
For a fascinating account of these struggles see Truesdell , The Rational 
Mechanics of Elastic or Flerible Bodies 1638-1788, L. Euleri Opera Omnia, 
Series II, Volume II, Part 2, Zurich, Fussli.? 

Consider a hanging cable with ends encased in rigid piers, as shown 
in Fig. 9.1. In the following, we let s denote distance along the cable 


1See the excellent, pioneering text by vanDyke, Perturbation Methods in Fluid Me- 
chanics, annotated edition, the Parabolic Press, 1975. 

2 As far as we know, the analogous problem for a cable with bending stiff ness-the 
full nonlinear version-was first set by one of us (jgs) on an examination in 1971 at 
the Technical University of Denmark. Shortly thereafter, a perturbation solution was 
published by A.M.A. van der Heijden, “On the Influence of the Bending Stiffness in 
Cable Analysis”, Proc. Kon. Ned. Ak. Wet. B76, 1973, pp. 217-229. 
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Figure 9.1: The geometry of and the forces on a hanging cable with stiffness. 


from its low point and w(s) its weight/length. Further, let H(s), V(s), and 
M(s) denote, respectively, the horizontal and vertical force and the moment 
exerted over the cross-section at s by the material to the right, as indicated 
in Fig. 9.1. Finally, with respect to a right-handed set of Cartesian axes 
Ozy, with O at the low point of the cable and the y-axis vertical, let z(s), 
y(s) and ¢(s) denote the Cartesian coordinates and tangent angle of a point 
P on the cable at a distance s. For the part of the cable between O and P, 
we have 


Horizontal force equilibrium: H(s) = H(0)= Ho. (9.1) 


Vertical force equilibrium: V(s) = | oUi: (9.2) 
Moment equilibrium: 
M(s) = M(0) — z(8)V (5) + y(s)H(s) + | “g(t)w(t)dt. (9.3) 
To complete our set of equations, we borrow from strength of materials 
the moment-curvature relation 
M(s) = EI(s)¢'(s) (9.4) 


where G = d¢/ds. Here, as with the drill string of Chapter 7, E is Young’s 
modulus and I is the moment of inertia of the cross-section. (Even though 
our problem is nonlinear, it may be shown that, so long as the cable is thin, 
a linear moment-curvature relation is quite accurate.) 
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We may reduce (9.1)—(9.4) to a second order DE for ¢ by differentiating 
(9.3), and, in the resulting equation, inserting (9.4) into the left side and 
(9.1) and (9.2) into the right side. Thus 


(r J eT (9.5) 
But 
sawdd, d single) . (96) 
Hence 
(EI¢’y + ( J l (9.7) 


We may restrict attention to the right half of the cable in which case (9.7) 
is to hold for 0 < s < L, where L is the length of the right half of the cable. 
The BC’s on G are then 


(0) = ¢(Z) =0. (9.8) 


Simplification and Nondimensionalization. Let us take E, J, and 
w to be constant and introduce the following dimensionless variable and 
parameters: 


2 . , k= wL/Ho , ê = EI/L? Hy. (9.9) 
Then (9.7) and (9.8) take the form 
620 + kz cos - sin G = O, 0< z<1, 600) = (1) O, (9.10) 


where a prime denotes differentiation with respect to z. In terms of L and 
$, E f 
128 J saih, yr f sin Gt) dt. (9.11) 
0 0 
The Catenary is the shape of a cable with no bending stiffness. We 
obtain the angle ® of the Catenary by setting « = 0 in the DE im (9.10): 


S Stan- (kz) . (9.12) 


As is obvious from the physics, ® fails to meet the BC at z = 1. Clearly 
our BVP has a singularity in the model. 


Exercise 9.1. Determine z and y for the Catenary from (9.11) 
and verify that y/L = cosh(z/L) — 1. 
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Boundary Layer Arguments. From experience, we know that if the 
cable is very thin (e? << 1), it must hang very nearly as a chain except 
near the ends. There ¢ must turn rapidly so that the cable can meet the 
pier horizontally. That is, ¢ must exhibit a boundary layer near z = 1. 
The dimensionless length 1 of the right half of the cable is one scale in our 
problem. The other is the dimensionless width of the boundary layer, say 
6. What is 6? We argue in rough terms as follows. From (9.12) we know 
that just outside the boundary layer, 6 W tan- i k = O(1). But as we 
move through the boundary layer, ¢ must decrease to zero to satisfy the 
BC. Thus ¢ undergoes a change of O(1) over a distance 6, which suggests 
that, if ¢ changes smoothly, then G in the boundary layer is O(5~"). Let 
us now integrate (9.10) across the boundary layer, from z = 1—6 to z = 1. 
This yields 


42%) e = f bine kes) dr Ole. (9.13) 


But outside the boundary layer the term ¢?¢’, which represents the effects 
of bending stress, is negligible. Thus, since we just concluded that ¢/(1) = 
O(ö5ᷣ—1), (9.13) tells us that ö = O(e). 


Boundary layer (or, more generally, order of magnitude) argu- 
ments are potent, heuristic guides for sorting out the various 
scales in a problem. Their justification is usually a posteriori, 
in the apparent consistency and reasonableness of the approxi- 
mations they suggest. 


The two-scale method may now be applied, systematically, to reduce 
our singular perturbation problem to a regular one. 
We introduce a fast (or boundary-layer) variable via 


1822 9.14) 


and assume that 


¢= o(¢, 2, €). (9.15) 
By the chain rule and (9.14), 


= —e7 1b + ¢:, G, e 2c = 26 1G + $z, (9.16) 
so that our BVP (9. 10) takes the form 


dee + kz cos - sin G — de: + epz = 0, OCC, 0<z<1, 
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G01, 0, €) = ¢(0, 1, €) = 0 " (9.17) 


We now assume that the introduction of the boundary-layer variable ¢ 
has made ¢ regular in e, i. e., that ¢ has the uniformly valid expansion 


5c. 2, 0 l C. Te C. ) 4 (9.18) 
Substituting (9.18) into (9.17), using the Taylor expansions 
9 i 93 3 0 
sin(g +e $ +- --) = sin G ted cos $ T (9.19) 
W 1 W a. fe 
cos(¢ +e t = cos G e & sin 6 (9.20) 


and equating to zero coefficients of like powers of e, we obtain the sequence 
of BVP's 


0 0 0 
ce +kz cos ġ sin ¢ = 0, b<f<e4, epg] 
$ (1,0) = (0,1) =0. (9.21a) 
1 0 o 1 0 
6 (Ez sin ¢ + cos 6) 6 -2¢, =0, O KCC, 21 
$ (e, 0) = (0,1) 0, (9.21b) 


etc. 


Exercise 9.2. Work out the equations coming after(9.21b). 


0 
A First Integral of the PDE in (9. 21a) is obtained by multiplying by ¢,. 
This enables us to rewrite the resulting equation as 


7050 + (Er sin $)c + (cos 9 =0, (9.22) 
which implies that 
1:9 0 0 
7 64 +kzsing + cosd=C(z),0<¢<6',0<2<1, (9.23) 
where C(z) is an unknown function of integration. 


The exact solution of (9.23) can be expressed in terms of elliptic inte- 
grals. However a remarkable simplification occurs (no elliptic integrals!) if 
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we first determine C(z). Imagine fixing a point P on the centerline of the 
cable and then letting the stiffness approach zero. Clearly, the slope of the 
centerline at P will approach that of the Catenary. That is, if we fix a 
value of z in (0,1) and let € — 0 then, from (9.12), ( œo and 

0 0 

6 (C. 2) ~ Bz), o&(¢, 2) ~ 0. (9.24) 


But (9.23) holds for all z in (0,1) and for all ¢ in (O, 1). Thus, in 
particular, it must hold in the limit as ¢ = e~! co, which implies that 


C(z) = kz sin (z) + cos z) (9.25) 
= sec & (z), 


where the last line comes from (9. 12). 
To further reduce (9.23) let 


6 (¢,2)=®(2) — . -) , 6 So. (9.26) 
Then, with the aid of a little trigonometry, we find that 
O¢ = —2sin(8/2)sect H(z) . (9.27) 


Exercise 9.3. Show the steps leading from (9.23) to (9.27). 


To solve (9.27), separate variables and integrate: 
J <=e(0/2)a(0/2) = —see4@(2) fac, (9.28) 


1. E., 
Intan(6/4) csec HHC) + D(z), 9.29) 
where D(z) is a function of integration. 
The first BC in (9. 21a) is satisfied automatically; the second, by (9.26), 
requires that 
60, 1) = ®(1) . (9.30) 
Hence, 
D(1) = Intan{;®(1)] (9.31) 


Setting 
D(z) = D(1) + E(z), (9.32) 
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solving (9.29) for 6, and substituting the resulting expression into (9.26), 
we obtain 


6 (C, z) = ®(z)—4tan-"{tan[®(1)/4] exp[—Csec/?6(z) + E(z)]}. (9.33) 
The function F(z) cannot be determined without considering the BVP 


(9.21b) for A However, since E(1) = O, we may ignore E(z) in (9.33) to 
0 


within the error made by approximating ¢ by ¢. To within this same error, 
we may also replace ®(z) in (9.33) by ®(1). 


Exercise 9.4. Justify these claims by showing that 
exp[—(sec/?6(z) + E(z)] = [1 + O] exp[—¢ sec? (i)]. 


With the aid of (9.12) and (9.14) and the approximations just men- 


tioned, we may express ¢ as follows: 


$~ tan-"(kz) — 
4 tan {tan( i tan k) exp a — wy] } . (9.34) 


A graph of the right side of (9.34) for two values of € and k = 0.6 is given 
in Fig. 9.2. 


Exercise 9.5 A designer is interested in the direct and bend- 
ing stresses in a cable, op and og. These are given by op = 
(H? + V?)/A and øg = (rE/L)d¢/dz, where A is the cross- 
sectional area and r is the radius of the cable. Taking Hp = 2000 
lbs., w = 5 lbs./ft.,r = 1 in., L = 200 ft., and E = 107 Ibs. /in. 2, 
compute the maximum values of op and og. 


A Problem from Cell Biology. Our modern example is one which 
certainly cannot be solved analytically and would be troublesome to solve 
numerically. The problem is to determine the effect of a protein on the 
diffusion of calcium through the walls of the intestines. We assume that the 

See Kretsinger, R.H., Mann, J.E., and Simmonds, J. G., “Evaluation of the Role of 


Intestinal Calcium Binding Protein in the Transcellular Diffusion of Calcium,” Proc. 5th 
Workshop on Vitamin D (A. W. Norman, Ed.) 1982, pp. 233-248. 
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Figure 9.2: Lowest approximation to the slope of the hanging cable for 
k = 0.6 and e = 0.1, 0.05. 


protein forms a compound with the calcium and that both the calcium and 
the calcium-protein compound are transported by diffusion. The calcium 
concentration is constant on the inside of the membrane and is pumped out 
by a Michalis-Menten pump at the outside membrane. The membrane is 
impervious to the protein: the protein cannot escape the cell. The following 
equations represent a mathematical formulation of the problem. 


dB I 
De r= i kon A-B + %% AB = O (9.35) 
ine — kon A-B + % AB=0 (9.36) 
AB — 
Dig σ + kon A-B — f AB =0 (9.37) 


dA d 


2 Ge =O ate=0ande=L. (9.39) 


In these equations, A is the concentration of protein [mol/cm*], B is the 
concentration of calcium, and AB is the concentration of the protein- 
calcium compound. The other parameters are identified in Table 9.1. 
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Table 9.1. Cell parameters 
Symbal 
3.16 x 10-*mol/cm” | calcium concentration at inside wall 
2.0 x 10-®cm?/sec diffusion constant for protein 
2.0 x 10-%cem?/sec diffusion constant for compound 
10-5em?/sec diffusion constant for calcium 


10h m/ mol sec rate of ſorming compound 

10?sec™! rate of decomposition of compound 
5 x 10-?mol/em?sec | max rate of Michaelis-Menten pump 
5 x 10 m thickness of membrane 


One final equation is required to specify the amount of protein in the 
membrane, namely 


Total Protein = H J 5 Adz + [ k Adr) (mol / em) (9.40) 


These equations are simplified by introducing dimensionless variables 
and parameters as follows: 


€=2/L, y= B/Bo, u=A/Bo, v =AB/Bo (9.41) 

À = VinaxL/Dp Bo, h = Km/Bo, ki BEV Da) en, k = ie 

o = Dg/Da, s= total protein / Bo. 445 

In (9.41) (9.43) we have taken Da = Dag. Equations (9.35) -(9. 39) now 
reduce to 

oy" — kjuy + kav =0 (9.44) 

u” — kjuy +kov =0 (9.45) 

vo” + kuy- kav =0 (9.46) 

ai a 

y(0)=1, -y (1)= 877 (9.47) 

1 (0) ) = O, u/(1)= (1) =0 (9.48) 

s= f udg + [ræ (9.49) 

0 0 


Equations (9. 44), (9. 45), (9.46) may be added to produce certain simpler 
equations, After some additions and integrations we arrive at 


ey — [1 + Ay — Cy'(0)ély - Cy? + Dis—7) + y (0) =0. (9.50) 
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In this equation A, C and D are constants whose values are related to the 
original parameters; 7, on the other hand, is a constant of integration from 
(9.44),(9.45), (9.46) which cannot be determined until the solution of (9.50) 
is known. Accompanying (9.50) are the conditions 


100) = 1, %% ), -y (1) =AvQ)/[(yG) + 4). (9.51) 


Equation the second equation in (9.51) is required to insure steady state 
conditions. From Table 9.1 we have 


6 = 8.0 * 1074, A = 0.632, C= 3.16, D = 0.200 . (9.52) 


If we set e = 0 in (9.50), we get the reduced equation 
— [1 + Ay — Cy (0X]Y - CY? + D(s— 7) + y (0) =0. (9.53) 


This equation is simply a quadratic in Y containing the two unknowns y 
and y/(0); we want the positive solution. To proceed, let 


Y(0) = L, (9.54) 


where L (which now stands for “left”) is the solution of (9.53) at € = 0. 
Near C = 0 we set 
y=Y + yee - (9.55) 


When (9.55) is substituted into (9.50), we find that 
eyar — [1 + Ay — Cy (oe + 2CY]ysz CY Te“ =0. (9.56) 


With 7 = €e~1/? as the boundary layer coordinate, we expand ygy in 
powers of 6/2 to obtain, to lowest order, 


n” =(1 + Ay + 2CL)n + Cn? (9.57) 
= P?Q?n + 505 R, 
here yas = N + em .. . . The solution of (9.57) which goes to zero as 


T goes to co is 
n= Presch'( Se Ue +a), (9.58) 


where we choose a so that 


1 — L=P*csch?a. (9.59) 
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Equation (9.59) insures that y(0) = 1. It now follows that to lowest order 
li. e., to within a relative error of O( el, 


7⁰⁰ K 202 (9.60) 
P Q(1— L)(1—L + P?)'/? 
= PA aa cosh a = I NE 


On the right end of the interval, we introduce the boundary layer coordinate 
r = (1 -C) e- 0/2 and again expand E in powers of / to obtain, to lowest 
order, 


= [1 + Aiy — Cy (0) + 2CR)\¢ - CC? (9.61) 
= PQ% 3050, 
where yg = —¢ + €1/2¢, + --- and Y(1) = R. We put the minus sign 


in the expansion of yg, to insure that Ç is positive. The solution which 
decays when € decreases is 


¢ = P’sec 212 O- 72(1— 6) + A). (9.62) 
The value of @ is obtained from 


y'(0) n 45 (== a FQ sech H sinh? . (9.63) 


Finally, the pump condition is satisfied if 
— ¥(0) = 9/(0) = A[R — ¢(1)]/[R - ¢(1) A (9.64) 


The unknowns, R, L, y enter (9.59), (9.60), (9.63), and (9.64) in some 
complicated, algebraic way. The easiest way to find the unknowns is in- 
versely. Proceed as follows. Choose a value for L (0 < L < 1). Use (9.51) 
with £ = 0 to express y in terms of this L. Use (9.60) to find y(0). Then £ 
can be determined from (9.63). Finally, (9.64) is satisfied by choosing À. In 
this procedure we must find R and check that R — ((I) is positive; if not 
then we need to choose a larger value for L. The final matched solution is 


y=Y +n-¢€ + O(e/?). (9.65) 


With the aid of this solution, we may plot curves like Figs. 9.3 and 9.4. 
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0.5 1 


Figure 9.3: Dimensionless concentration of calcium vs. distance across the 
cell. 


0 8 10 15 20 


Figure 9.4: Dimensionless flux of calcium vs. pump rate À for various 
values of the dimensionless cell protein s. 
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The above problem, though complicated, illustrates the power of pertur- 
bation theory with matched expansions to solve nonlinear boundary value 
problems. The two-scale method of Chapter 5 could also have been used to 
obtain an approximate solution to the Calcium diffusion problem and we 
urge to reader to do so. 
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Appendix A 


Roots of Te(z) and Tp(z) 


The lemma on page 33 means: given any root 25 of To(z) of multiplicity 
yk (Yı + 72 +--+ = n), and given any constant p > 0 no greater than 
half the distance from z to the nearest distinct root of Ty(z) (call it zm), 
there exists a value of € such that there are precisely 7 roots of T,(z) in 
a disk of radius p centered at 25. See Fig. A. 1. A student who has taken 
a first course in complex variable theory should have no trouble with the 
following proof. 

Given a function f(z) that is continuous and non-vanishing on, and 
analytic within, a simple, rectifiable, closed curve T, the number of zeros 
of f within T is, by Rouché’s Theorem, 


N(f) = = 5 700 dz, (A.) 


where F (z) = df/dz. The idea is to show that if we take T to be a circle 
of radius p centered at 25, then 


kans To(z)[To(z) + Ee(z)] d (A.2) 


INT) — NTO) 25 


will be less than 1 if e is sufficiently small. We do this by obtaining a lower 
bound on the magnitude of the denominator in (A.2) that does not vanish 
with € and an upper bound on the magnitude of the numerator that does. 
The key is to note that the hypothesis lim E.(z) = 0 as € 0 implies two 
things: First, if e is sufficiently small, 


|E<(z)| <(1/2)p” , z eT, (A.3) 
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Figure A.1: Some roots of T. (z) and T(z). 


since 2 is bounded on T. Second, as E,(z) is of the form 
E.{z) = Cm(e)z™ . . + Ci (e) z + Cole), (A. 4) 
then, with C(e) = max|(C;(€)|, j = 1, 2, m and r = zz, 


Ee( z) < |Cm(e)llr + f +---+ |Co(€)| 
< mC(|\r + f +1), 2 ST. (A.5) 


By the same reasoning, 
IEN) < m(m - od = +1), zer. (AS) 
To obtain a lower bound on T(z) note from its factored form 
m)l = lz alls - zll- zal (A.7) 


that 
ITo(z)|> P , zer. (A.8) 


To obtain a n upper bound on To(z) and 78 (z), let R = max{|z;—zz|, 1}, j = 
1,2, ... n. Then (A.7) implies that 


To(z) < (R + p), z ET. (A.9) 
Furthermore, as 


Tels) = ( -aa + C aE hji — 0) 
+e- + (z — 21) (2 — 28-1), (A. 10) 
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and as there are n terms each containing n — 1 factors, 
IT| S NR + A! , 2 ST. (A.11) 
Thus in (A.2) 
ToC z) ToC) + Ee(z)| > NTC) - lEz) (1/2) (A.12) 
and 


|To(2)Ei(z) — E,(z)T6(2)| 
< [T(E + IET) (A.13) 
< C(e)[m(m — „eh e + 1)(R + p)” 

+ mn(|r + pl” + DR + Ga-]. 


Thus 


|N(T.) - N(D)| 
< 2p'-?"C(e)m(R + p)" im — 1) 
x(n f + 1 R + p) + n(lr + ol}. (A.14) 


But lim Ee(z) = 0 as e — 0 implies that lim C(e) = 0. Thus, given any 
p > 0, the right side of (A.14) is less than 1 if e if sufficiently small. Q. E. D. 


Appendix B 


Proof that 
Ry, = OT] 


To prove that (4.35) is implied by the IVP (4.32), we first replace v in 
(4.32) by the right side of (4.34). Then, noting that A has the form (4.33) 
and that 20, „z satisfy the sequence of equations (4.36), we obtain an 


IVP for Ry +; of the form 
X Rygi RAI = Af (z0, zv, Nei), RA) = RVO) = 0. (B. I) 


By the method of variation of parameters, we convert (B.1) into the non- 
linear integral equation 


Ni (T, B) = BA" fi sin- (T — 7) f(z0,--+,Rw4i(7,8))dr. (B. 2) 


The potential energy associated with (4.32) is 1/221 (1/12) BJ. Hence 
(4.32) has periodic solutions if Ø < 6. By the choice of À this period is 
2x. But 20, „ z are, by construction, 27-periodic functions. Hence, by 
(4.34), Ri must be 27-periodic. Thus, in (B.2), we may, without loss of 
generality, assume that 0 < T < 27. 

The general method of proof may be inferred from that for N = 0, which 
we choose to minimize algebraic complexity. Thus, with A? = 1 + By(f), 
u(B) = O(1) and 20 = cos T, (B.2) reduces to 


Ri (T. 8) = BA? [ " sin A (T- 7){a(8)cosT + [cost + RIC, ) Idr. 
(B.3) 
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Taking absolute values of both sides of (B.3) and noting that all trigono- 
metric functions are bounded in absolute value by 1, we infer the inequality 


R = |Rı(T,2)| < BA? | “tl + (14+ R)’]dr (B.4) 
$ T 
SBR (14 hear, 


where k = 1 + |a|. Now (B. 4) implies that R(T) < S(T), 0 < T < 2r, 
where 


S HN, S(0) = 0. (B.5) 


The IVP (B.5) may be solved readily, in closed-form, by separation of 
variables. This yields 


R(T) S ~1 + (1-21 17) 72. (B.6) 


But k, A, and T are O(1). Thus, by (B.6) and periodicity, R(T) = 
|Ri(T, B)| = O(2), for all T. Q.E.D. 


Appendix C 


Approximate Evaluation 
of Integrals 


Certain integrals can be evaluated using techniques related to the ones we 
have used for differential equations. Such integrals usually arise when using 
integral transform methods to solve ordinary or partial differential equa- 
tions. We will demonstrate several techniques without actually deriving 
the integral itself. We will consider Laplace integrals first, then Fourier 
integrals, and finally integrals that can be approximated using integration 
by parts. 
Laplace integrals have the form 


I(s) = y f(t)e~ dt, (C.1) 


where f(t) meets the conditions required for I to exist. Thinking about the 
exponential function for s >> 1, we realize that the value of e“ is nearly 
zero unless t is close to zero. We would expect that the value of the integral 
would get smaller as s gets larger. Also the behavior of f(t) near t = 0 will 
be of primary importance because e~** is nearly zero for t larger than s~?. 


See Fig C.1. 
To see how I(s) behaves, we assume that f(t) may be expressed as 


f(t) = ao Tait +--+ + ant” + O(t"**) for O < t ST. (C.2) 
Noting that f° = f + [© and inserting (C. 2) into (C. I), we obtain 
I(s) = J : e [ao + ait +--+ ant” + O(t"*)]dt + J x e (t) dt 
0 T 
= J(s)+ K(s). (C.3) 
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We shall show that K(s), the second integral in (C.3), is transcendentally 
small and then use the first integral, J(s), as our approximation for I(s). 
One requirement for (C.1) to exist is that there exist numbers M and b 


such that 
If (t)| < Me™ for T <t. (C.4) 


o0 o⁰ M 
pi 0 ef e" Me” dt = „ , 
T ~ s—b 


which shows that the second integral in (C.3) is transcendentally small. 
For J (s), the first integral in (C.3), we need to recall that O(¢"*") in (C.2) 
means that there exists some positive constant A such that 


[O(t S AH for 0 <t <T. 


Therefore 


— 
<A ji eta dt. (0.5) 
0 


J(s) = D ay p en *tt* dt 


0 


We notice that all integrals in the inequality (C. 5) may have their limits 
extended to oo by introducing terms that are transcendentally small. This 
is seen by noting that 


K oo oo k! 
/ etats | e dt = e ft Err 
0 y 0 8 


Moreover, - = 
[ e tk dt = oe] (7 ＋ T) dr, 
T 0 


fit) 


-St 


Figure C.1: The integrand of (C.1) showing f(t) and the rapid decay of the 
exponential. 
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as may be seen by setting t = T +T. The right side of this last equation 
goes to zero exponentially fast as s — oo and is therefor transcendentally 
small. Putting the above pieces together, we may write 


œ ‘ nl * 
aA 24 S44 So (C8) 


The proof above may be extended to integrals of the form 
een, a>- 
0 
where f(t) satisfies (C. 2) and (C. 4). Using the gamma- function 


oo 
T(z) = f ett. Id, 
0 
and a proof similar to the one above, we may establish 


Watson’s Lemma: Under appropriate conditions on f(t) for the integral 
to exist, 


= T(o+1) , al(a+2 
i emt F(t)dt aaa = 3 = = 


a,T'(a+n+1) 
2 satn+l , 


as s — œo and the a;’s are the coefficients in (C.2). 


a>-l, 


Example: Find the expansion for 


= [ nie a 
as §— O. 


Solution: Notice that the integrand is in the form of the one in Watson’s 
Lemma with a = 1/2. Also 


fi) = = 74 . 1-47). 


at 


(Although the series diverges for 1 < t, f(t) has the form of (C.2) for any 
fixed T > 0.) Using Watson’s Lemma, we have 


T(3/2) 17(67%2 
I(s) = CP) PUM) o (gr) 1 
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Laplace’s method is based on the above expansion for Laplace inte- 
grals but is used for a broader class of integrals of the form 


I(s) = / Faje dt, (C.7) 


where we suppose that g(t) has a minimum at to € [a,b]. From what has 
been said above, we expect that the major contribution to I(s) will come 
from a small interval containing to. In this interval, we replace f and g by 
suitable approximations and use Watson’s Lemma to secure an approxima- 
tion for I(s). If to is at an end point of the interval of integration and if 
f'(to) 2 0, then the integral is approximated by a straight forward exten- 
sion of Watson’s Lemma. Otherwise, we will assume that g(t) has a Taylor 
series of the form 


a 
t 
a(t) = gto) + CEL- ta)? +--+ (C8) 
with g”(to) > 0. If g has a minimum whose Taylor series does not have 
this form, the procedure will have to be modified. Since the behavior of g 
is important only in a small neighborhood of to, we replace (C.7) with 
In~h= [ Ff (them toto + g-t (C9) 
a 
By letting g’(to)(t — to)?/2 = u, we can write 


b 
Ty = e 0 /2/g" (to) f | f(to + /2/9"(to)u)e-** du (C.10) 
where a* = \/g"(to)/2(a - to) and b* = yg” (to)/2)(b - to). Since a* < 0 


and ö“ > 0. We incur only transcendentally small errors if we take the 
upper and lower limits in (C.10) to be oo and —oo respectively. By letting 
u = r for u > 0 and u = M for u < 0, we can reduce (C.10) to 


= 1/2e~#9(t0) /2]9"(to) „17 f(to+ [(2] 9" (to)x)e~** dr r 
+ o f(to - Vo) a e (C. 11) 
Each of the integrals in (C. 11) may be approximated by using Watson’s 


Lemma. In a specific problem, we usually deal with the functions at hand 
rather than use the formalism above. 
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Ezample: Find an expansion for 


I(s) = J e. 


Solution: The integral is not in the precise form of (C. 7), but s cost = 
—s(— cos t) fixes the defect. We then note that — cos t has its minimum 
value at t = 0. We make the change of variable cost = 1 — u? (we are 
making an exact change of variable, not approximating cost), and obtain 


a ae 
I(s) f em rr Ta aii 
where a = — y1 —cos(—1) and b* = \/1 —cos(1.2). If we change a* to 


—oo and b* to oo, the error is transcendentally small. The integrand is 
now an even function of u and the limits of integration are symmetric with 
respect to the origin. Therefore, we may double the value of the integral 
and change the limits of integration to obtain 

2du 
uz 
Making the interval of integration contain only positive numbers makes it 
simpler to make the substitution u? = z and yields 


I(s) c- [A r 97 (C.12) 


We can readily apply Watson’s Lemma to (C.11) when we use the binomial 
expansion to write 


f(z) = 2-20 1 


I(s) = 2e’ ow 
0 


= hraa (3) ()) 
aan ae * J. (C.13) 


Thus with a = —1/2, we find 


T(1/2) T(3%½) 31652) 
I(s) ~ ~ Vie | 8177 2 re —＋ 


Exercise C. 1. Find an asymptotic expansion for s > 1: 
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oo 
(a) T(s+ 1) = e~*t*dt 
Hint: Let t a su. When you have completed the prob- 
lem, you will have derived a version of Stirling’s for- 
mula. 


(b) Ie) = f id e. cot dt. 


1.2 
© a) = f ena, 
-1 
Note that the solution process for part (b) will be quite 
different from that of part (c). How does the fact that 
the integrand of part (c) is odd affect the solution? 
= dt 
0 


(d) I(s) sf * 


Fourier Integrals, another import ant class of integrals which may have 
an asymptotic expansion for s = w >> 1 are of the form 


fis) 2 i f 7000 ar ae, \ ae (C.14) 


The Riemann-Lebesgue Lemmastates that lim. J(w) = 0 provided that 
f(z) and f?(z) are integrable on (a,b). The lemma is proved in many 
advanced calculus and analysis books. (e.g., Fulks, W., Advanced Calculus, 
3rd ed., p.651, Wiley, 1978.) If f has several continuous derivatives, we 
may use the Riemann-Lebesgue Lemma in conjunction with integration by 
parts to find an asymptotic expansion for J(w). Choosing the cos in (C.14), 
we have 


Hr = [ " f(2)sin(we)de (015) 


w 


[In passing, we note that (C.15) could be used to prove the Riemann- 
Lebesgue Lemma if the integral is bounded.] In (C.15), we use the Riemann- 
Lebesgue Lemma to conclude the first term on the right is larger than the 
second for w >> 1 since the integral goes to zero and is multiplied by w~?. 
The process of integration by parts may be continued until some derivative 
of f becomes so badly behaved that the Riemann-Lebesgue Lemma no 
longer holds. In this way, we obtain the expansion of f in powers of w~?. 
Example: Find an asymptotic expansion for 


I(w) = ; . 
3 
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Solution: Since (1 + 22) 1 has derivatives of all orders, we may integrate 
by parts as many times as we wish to obtain a series with any number of 
terms. The first integration produces 


(1/26) sin(5w) — (1/10) sin(3w) | 1 2z sin( or) 2z sin(wz) 4 


I(w) = * G42) 


If we carry the process on for one more step, we find 


Iw) ——— (1/10) sin(3w) 
+ (8/50) cos(3w) — (6/338) cools) 1 /> (2-62?) cos(wz) 4, 


w? Ta 3 (l + z?)8 


The terms in the equation above are in increasing powers of w~! and the 
last term is the smallest since the integral itself goes to zero as w — oo. 

Though we will not prove the Riemann-Lebesgue Lemma, we do wish 
to comment on why the theorem works: when w is very large, the cos(wz) 
term goes through a complete period in a very small portion of the interval 
of integration. Moreover, ie = cos(wx)dz = 0 for any a in the interval 
of integration. Thus, if w is very large, the function does not have a chance 
to change much while cos(wz) is producing canceling areas. The major 
contribution to the integral comes from the end points of the interval, as 
is apparent from (C.15). This idea of rapid oscillation is the key idea in a 
method used in conjunction with complex integration called the method of 
stationary phase. 


Exercise C.2. Find asymptotic expansions for 


(a) / ne sin(wz)dz. 


-4 
Find the exact value of the integral and compare the 
value to the asymptotic value for w = 100 and 1000. 


J dz 
b) J cor) ==. 
With an appropriate change of variable, this integral 
can be expressed in terms of the Fresnel integral. See 
Abramowitz and Stegun. 


Integration by parts is the last method we will mention. This method 
is easy to use; indeed, we havealready used it for the trigonometric integrals 
above and for the Stieltjes’ integral in Chapter 1, (1.69). However, it is 
hard to give a general theory for the method and it is hard to describe 
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those situations in which the method will either work or fail. Since we have 
given examples previously, we will simply set the following integrals for you 
to do. 


Exercise C.3. Find asymptotic expansions for the following 
integrals. Assume that À is large in all cases. You must show 
that the integral which is left is smaller than the other parts of 


the expression. 


oa e 
(a) 10. f = 
* 
(b) An integral related to the incomplete T-Function is 


given by 
yla, ) = £ e de: 
0 


(e) The complimentary error function 


erfc(A) = 7 e`" dt 
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